#*-------------------------------------------------------------------
* 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: Paula B. Staudt
* $Id: tray.mso 522 2008-05-21 23:21:12Z arge $
*--------------------------------------------------------------------*#

using "streams";

Model tray

ATTRIBUTES
	Pallete 	= false;
	Icon 		= "icon/Tray"; 
	Brief 		= "Complete model of a column tray.";
	Info 		=
"== Assumptions ==
* both phases (liquid and vapour) exists all the time;
* thermodymanic equilibrium with Murphree plate efficiency;
* no entrainment of liquid or vapour phase;
* no weeping;
* the dymanics in the downcomer are neglected.

== Options ==
You can choose the equation for the liquid outlet flow and the vapour
inlet flow calculation through the VapourFlowModel and LiquidFlowModel
switchers.

== References ==
* ELGUE, S.; PRAT, L.; CABASSUD, M.; LANN, J. L.; CéZERAC, J. Dynamic models for start-up operations of batch distillation columns with experimental validation. Computers and Chemical Engineering, v. 28, p. 2735-2747, 2004.
* FEEHERY, W. F. Dynamic Optimization with Path Constraints. Tese (Doutorado) - Massachusetts Institute of Technology, June 1998.
* KLINGBERG, A. Modeling and Optimization of Batch Distillation. Dissertação (Mestrado) - Department of Automatic Control, Lund Institute of Technology, Lund, Sweden, fev. 2000.
* OLSEN, I.; ENDRESTOL, G. O.; SIRA, T. A rigorous and efficient distillation column model for engineering and training simulators. Computers and Chemical Engineering,v. 21, n. Suppl, p. S193-S198, 1997.
* REEPMEYER, F.; REPKE, J.-U.; WOZNY, G. Analysis of the start-up process for reactive distillation. Chemical Engineering Technology, v. 26, p. 81-86, 2003.
* ROFFEL, B.; BETLEM, B.; RUIJTER, J. de. First principles dynamic modeling and multivariable control of a cryogenic distillation column process. Computers and Chemical Engineering, v. 24, p. 111-123, 2000.
* WANG, L.; LI, P.; WOZNY, G.; WANG, S. A start-up model for simulation of batch distillation starting from a cold state. Computers and Chemical Engineering, v. 27, p.1485-1497, 2003.
";	
	
PARAMETERS
outer PP 		as Plugin		(Brief = "External Physical Properties", Type="PP");
outer NComp 	as Integer 		(Brief="Number of components");

	Mw(NComp)  							as molweight 		(Brief="Component Mol Weight",Hidden=true);
	Gconst 								as acceleration 	(Brief="Gravity Acceleration",Default=9.81,Hidden=true);
	zero_flow 							as flow_mol 		(Brief = "Stream Flow closed",Default = 0, Hidden=true);
	low_flow  							as flow_mol 		(Brief = "Low stream Flow",Default = 1E-6, Hidden=true);
	Pi									as constant 		(Brief="Pi Number",Default=3.14159265, Symbol = "\pi",Hidden=true);
	
	TrayDiameter_			as length 		(Brief="Tray Diameter",Default=1.600);
	TraySpacing_			as length 		(Brief="Tray Spacing",Default=0.600);
	Fraction_HoleArea_ 		as fraction 	(Brief="Fraction of the active area that is occupied by the holes with respect to the total tray area",Default=0.10);
	Fraction_DowncomerArea_ as fraction 	(Brief="Fraction of the downcomer area with respect to the total tray area",Default=0.20);
	WeirLength_				as length 		(Brief="Weir length", Default = 1);
	WeirHeight_ 			as length 		(Brief="Weir height", Default= 0.05);
	TrayLiquidPasses_ 		as positive 	(Brief="Number of liquid passes in the tray", Lower = 1,Default=1);
	HeatSupply_ 			as heat_rate 	(Brief="Rate of heat supply",Default = 0);
	AerationFraction_ 		as Real 		(Brief="Aeration fraction", Default = 1);
	DryPdropCoeff_ 			as Real 		(Brief="Dry pressure drop coefficient", Default= 0.60);
	MurphreeEff_ 			as Real 		(Brief="Murphree efficiency for All Trays",Lower=0.01,Upper=1);

	PlateArea_ 				as area 		(Brief="Plate area = Atray - Adowncomer",Protected=true);
	TrayVolume_				as volume		(Brief="Total Volume of the tray",Protected=true);
	HolesArea_ 				as area 		(Brief="Total holes area",Protected=true);
	
	FeeheryCoeff	as Real 	(Brief="Feeherys correlation coefficient", Unit='1/m^4', Default=1,Hidden=true);
	ElgueCoeff		as Real 	(Brief="Elgues correlation coefficient", Unit='kg/m/mol^2', Default=1,Hidden=true);
	OlsenCoeff		as Real 	(Brief="Olsens correlation coefficient", Default=1,Hidden=true);
	
	VapourFlow 	as Switcher	(Brief="Flag for Vapour Flow condition",Valid = ["on", "off"], Default = "off",Hidden=true);
	LiquidFlow 	as Switcher	(Brief="Flag for Liquid Flow condition",Valid = ["on", "off"], Default = "off",Hidden=true);
	
VARIABLES

	Inlet 				as stream 			(Brief="Feed stream", Hidden=true, PosX=0, PosY=0.4932, Symbol="_{in}");
	LiquidSideStream	as liquid_stream 	(Brief="liquid Sidestream", Hidden=true, Symbol="_{outL}");
	VapourSideStream 	as vapour_stream 	(Brief="vapour Sidestream", Hidden=true, Symbol="_{outV}");

in	InletLiquid 	as stream 			(Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
in	InletVapour 	as stream 			(Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
out	OutletLiquid 	as liquid_stream 	(Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
out	OutletVapour 	as vapour_stream 	(Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");

	LFlowModel 	as positive	(Brief="Flag for Liquid Flow Model",Lower = 0, Default = 1 , Hidden=true);
	VFlowModel 	as positive	(Brief="Flag for Vapour Flow Model",Lower = 0, Default = 1 , Hidden=true);

	M(NComp) 			as mol 			(Brief="Molar Holdup in the tray");
	ML 					as mol 			(Brief="Molar liquid holdup");
	MV 					as mol 			(Brief="Molar vapour holdup");
	E 					as energy 		(Brief="Total Energy Holdup on tray");
	vL 					as volume_mol 	(Brief="Liquid Molar Volume");
	vV 					as volume_mol 	(Brief="Vapour Molar volume");
	Level 				as length 		(Brief="Height of clear liquid on plate");
	yideal(NComp) 		as fraction;
	rhoL 				as dens_mass 	(Brief="Mass Density of liquid phase");
	rhoV 				as dens_mass	(Brief="Mass Density of vapour phase");

SET

	Mw = PP.MolecularWeight();
	zero_flow = 0 * 'kmol/h';
	low_flow = 1E-6 * 'kmol/h';

	PlateArea_ = 0.25*Pi*(TrayDiameter_^2)*(1-Fraction_DowncomerArea_);
	TrayVolume_ = 0.25*Pi*(TrayDiameter_^2)*TraySpacing_;
	HolesArea_ = 0.25*Pi*(TrayDiameter_^2)*Fraction_HoleArea_;

EQUATIONS

# LiquidFlow and VapourFlow equations need to be linerized to avoid indetermination !
switch LiquidFlow

	case "on":
			
			if LFlowModel equal 1 then
				"Francis Equation"
				OutletLiquid.F*vL = 1.84*'1/s'*WeirLength_*((Level-(AerationFraction_*WeirHeight_))/(AerationFraction_))^2;
			
			else if	LFlowModel equal 2 then 
				"Wang_Fl"
				OutletLiquid.F*vL = 1.84*'m^0.5/s'*WeirLength_*((Level-(AerationFraction_*WeirHeight_))/(AerationFraction_))^1.5;
			
			else if	LFlowModel equal 3 then
				"Olsen"
				OutletLiquid.F / 'mol/s'= WeirLength_*TrayLiquidPasses_*rhoL/sum(Mw*OutletVapour.z)/(0.665*OlsenCoeff)^1.5 * ((ML*sum(Mw*OutletLiquid.z)/rhoL/PlateArea_)-WeirHeight_)^1.5 * 'm^0.5/mol';
			
			else if	LFlowModel equal 4 then  
				"Feehery_Fl"
				OutletLiquid.F = WeirLength_*rhoL/sum(Mw*OutletLiquid.z) * ((Level-WeirHeight_)/750/'mm')^1.5 * 'm^2/s';
			
			else   
				"Roffel_Fl"
				OutletLiquid.F = 2/3*rhoL/sum(Mw*OutletLiquid.z)*WeirLength_*(ML*sum(Mw*OutletLiquid.z)/(PlateArea_*1.3)/rhoL)^1.5*sqrt(2*Gconst/
							(2*(1 - 0.3593/'Pa^0.0888545'*abs(OutletVapour.F*sum(Mw*OutletVapour.z)/(PlateArea_*1.3)/sqrt(rhoV))^0.177709)-1)); #/'(kg/m)^0.0888545/s^0.177709';
end
end
end
end
		
		when Level < (AerationFraction_ *WeirHeight_) switchto "off";
		
		case "off":
		
		"Low level"
		OutletLiquid.F = zero_flow;
		
		when Level > (AerationFraction_ * WeirHeight_) switchto "on";
		
end
	
switch VapourFlow
	
	case "on":
			
			if VFlowModel equal 1 then
				"Reepmeyer"
				InletVapour.F*vV = sqrt((InletVapour.P - OutletVapour.P)/(rhoV*DryPdropCoeff_))*HolesArea_;
			
			else if VFlowModel equal 2 then
				"Feehery_Fv"
				InletVapour.F = rhoV/PlateArea_/FeeheryCoeff/sum(Mw*OutletVapour.z) * sqrt(((InletVapour.P - OutletVapour.P)-(rhoV*Gconst*ML*vL/PlateArea_))/rhoV);
			
			else if VFlowModel equal 3 then
				"Roffel_Fv"
				InletVapour.F^1.08 * 0.0013 * 'kg/m/mol^1.08/s^0.92*1e5' = (InletVapour.P - OutletVapour.P)*1e5 - (AerationFraction_*sum(M*Mw)/(PlateArea_*1.3)*Gconst*1e5) * (rhoV*HolesArea_/sum(Mw*OutletVapour.z))^1.08 * 'm^1.08/mol^1.08';
			
			else if VFlowModel equal 4  then
				"Klingberg"
				InletVapour.F * vV = PlateArea_ * sqrt(((InletVapour.P - OutletVapour.P)-rhoL*Gconst*Level)/rhoV);
			
			else if VFlowModel equal 5 then
				"Wang_Fv"
				InletVapour.F * vV = PlateArea_ * sqrt(((InletVapour.P - OutletVapour.P)-rhoL*Gconst*Level)/rhoV*DryPdropCoeff_);
				
			else 
				"Elgue"
				InletVapour.F  = sqrt((InletVapour.P - OutletVapour.P)/ElgueCoeff);
end
end
end
end
end
	
	when InletVapour.F < low_flow switchto "off";

	case "off":
		InletVapour.F = zero_flow;
		
	when InletVapour.P > OutletVapour.P switchto "on";

end

"Murphree Efficiency"
	OutletVapour.z =  MurphreeEff_ * (yideal - InletVapour.z) + InletVapour.z;
	
"Energy Balance"
	diff(E) = ( Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h- OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h
	- VapourSideStream.F*VapourSideStream.h - LiquidSideStream.F*LiquidSideStream.h + HeatSupply_ );

"Energy Holdup"
	E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletLiquid.P*TrayVolume_;

"Geometry Constraint"
	TrayVolume_ = ML* vL + MV*vV;

"Level of clear liquid over the weir"
	Level = ML*vL/PlateArea_;

"Component Molar Balance"
	diff(M)=Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z-
	LiquidSideStream.F*LiquidSideStream.z-VapourSideStream.F*VapourSideStream.z;

"Molar Holdup"
	M = ML*OutletLiquid.z + MV*OutletVapour.z;

"Mol fraction normalisation"
	sum(OutletLiquid.z)= 1.0;

"Mol fraction constraint"
	sum(OutletLiquid.z)= sum(OutletVapour.z);

"Liquid Volume"
	vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Vapour Volume"
	vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z);

"Liquid Density"
	rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);
	
"Vapour Density"
	rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z);

"Chemical Equilibrium"
	PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, yideal)*yideal;

"Thermal Equilibrium"
	OutletVapour.T = OutletLiquid.T;

"Mechanical Equilibrium"
	OutletVapour.P = OutletLiquid.P;
	
"Thermal Equilibrium Vapour Side Stream"
	OutletVapour.T = VapourSideStream.T;

"Thermal Equilibrium Liquid Side Stream"
	OutletLiquid.T = LiquidSideStream.T;

"Mechanical Equilibrium Vapour Side Stream"
	OutletVapour.P= VapourSideStream.P;

"Mechanical Equilibrium Liquid Side Stream"
	OutletLiquid.P = LiquidSideStream.P;

"Composition Liquid Side Stream"
	OutletLiquid.z= LiquidSideStream.z;
	
"Composition Vapour Side Stream"
	OutletVapour.z= VapourSideStream.z;

end

#*-------------------------------------------------------------------
* Model of a tray with reaction
*-------------------------------------------------------------------*#
Model trayReac
	ATTRIBUTES
	Pallete 	= false;
	Icon 		= "icon/Tray"; 
	Brief 		= "Model of a tray with reaction.";
	Info 		=
"== Assumptions ==
* both phases (liquid and vapour) exists all the time;
* thermodymanic equilibrium with Murphree plate efficiency;
* no entrainment of liquid or vapour phase;
* no weeping;
* the dymanics in the downcomer are neglected.
	
== Specify ==
* the Feed stream;
* the Liquid inlet stream;
* the Vapour inlet stream;
* the Vapour outlet flow (OutletVapour.F);
* the reaction related variables.
	
== Initial ==
* the plate temperature (OutletLiquid.T)
* the liquid height (Level) OR the liquid flow OutletLiquid.F
* (NoComps - 1) OutletLiquid compositions
";

PARAMETERS

	outer PP 			as Plugin(Type="PP");
	outer NComp 	as Integer;
	
VARIABLES

	Inlet 						as stream 				(Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
	LiquidSideStream	as liquid_stream 	(Brief="liquid Sidestream", Hidden=true, Symbol="_{outL}");
	VapourSideStream as vapour_stream 	(Brief="vapour Sidestream", Hidden=true, Symbol="_{outV}");
	

in		InletLiquid 		as 	stream 				(Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
in		InletVapour 		as 	stream 				(Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
out	OutletLiquid 	as 	liquid_stream 		(Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
out	OutletVapour 	as 	vapour_stream 	(Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");

	yideal(NComp) 	as fraction;

	M(NComp) 		as mol 				(Brief="Molar Holdup in the tray");
	ML 						as mol 				(Brief="Molar liquid holdup");
	MV 						as mol 				(Brief="Molar vapour holdup");
	E 						as energy 			(Brief="Total Energy Holdup on tray");
	vL 						as volume_mol 	(Brief="Liquid Molar Volume");
	vV 						as volume_mol 	(Brief="Vapour Molar volume");
	Level 					as length 			(Brief="Height of clear liquid on plate");
	Vol 						as volume;
	
	rhoL 			as dens_mass;
	rhoV 			as dens_mass;
	r3 				as reaction_mol 	(Brief = "Reaction resulting ethyl acetate", DisplayUnit = 'mol/l/s');
	C(NComp) 	as conc_mol 		(Brief = "Molar concentration", Lower = -1); 
	
EQUATIONS

"Molar Concentration"
	OutletLiquid.z = vL * C;
	
"Reaction"
	r3 = exp(-7150*'K'/OutletLiquid.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4))*'l/mol/s';
	
"Molar Holdup"
	M = ML*OutletLiquid.z + MV*OutletVapour.z;

"Thermal Equilibrium Vapour Side Stream"
	OutletVapour.T = VapourSideStream.T;

"Thermal Equilibrium Liquid Side Stream"
	OutletLiquid.T = LiquidSideStream.T;

"Mechanical Equilibrium Vapour Side Stream"
	OutletVapour.P= VapourSideStream.P;

"Mechanical Equilibrium Liquid Side Stream"
	OutletLiquid.P = LiquidSideStream.P;

"Composition Liquid Side Stream"
	OutletLiquid.z= LiquidSideStream.z;
	
"Composition Vapour Side Stream"
	OutletVapour.z= VapourSideStream.z;

"Mol fraction normalisation"
	sum(OutletLiquid.z)= 1.0;

"Liquid Volume"
	vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Vapour Volume"
	vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z);

"Thermal Equilibrium"
	OutletVapour.T = OutletLiquid.T;

"Mechanical Equilibrium"
	OutletVapour.P = OutletLiquid.P;

	Vol = ML*vL;
	
"Liquid Density"
	rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Vapour Density"
	rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z);

"Chemical Equilibrium"
	PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = 	PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, yideal)*yideal;

	sum(OutletLiquid.z)= sum(OutletVapour.z);

end

#*-------------------------------------
* Model of a packed column stage
-------------------------------------*#
Model packedStage

ATTRIBUTES
	Pallete 	= false;
	Brief 		= "Complete model of a packed column stage.";
	Info 		=
"== Specify ==
* the Feed stream
* the Liquid inlet stream
* the Vapour inlet stream
* the stage pressure drop (deltaP)
	
== Initial ==
* the plate temperature (OutletLiquid.T)
* the liquid molar holdup ML
* (NoComps - 1) OutletLiquid compositions
";

PARAMETERS

outer PP 	as Plugin	(Brief = "External Physical Properties", Type="PP");
outer NComp as Integer	(Brief = "Number Of Components");

	LiquidResistanceCoeff 	as positive 	(Brief="Resistance coefficient on the liquid load", Default=1,Hidden=true);
	AreaPerPackingVolume 	as Real 		(Brief="surface area per packing volume", Unit='m^2/m^3',Hidden=true);
	ColumnInternalDiameter	as length 		(Brief="Column diameter",Hidden=true);
	PackingVoidFraction		as Real 		(Brief="Void fraction of packing, (m^3 void space/m^3 packed bed)",Hidden=true);
	HeightOfPacking 		as length 		(Brief="Height of packing",Hidden=true);
	Number_Stages 			as Integer		(Brief="Number of Stages", Default=3,Hidden=true);
	HeatOnStage 			as heat_rate 	(Brief="Rate of heat supply",Hidden=true);

	HETP 			as length 		(Brief="The Height Equivalent to a Theoretical Plate",Hidden=true);
	ColumnArea 		as area			(Brief="Column Sectional Cross Area",Hidden=true);
	V 				as volume		(Brief="Total Volume of the tray",Hidden=true);
	Pi				as constant 	(Brief="Pi Number",Default=3.14159265, Symbol = "\pi",Hidden=true);
	Gconst 			as acceleration (Brief="Gravity Acceleration",Default=9.81,Hidden=true);
	
	low_flow  		as flow_mol 	(Brief ="Low Flow",Default = 1E-6, Hidden=true);
	low_pressure  	as pressure 	(Brief ="Low Pressure",Default = 1E-6, Hidden=true);
	zero_flow 		as flow_mol 	(Brief ="No Flow",Default = 0, Hidden=true);

	Mw(NComp)  	as molweight 	(Brief = "Component Mol Weight",Hidden=true);
	VapourFlow  as Switcher		(Brief = "Vapour Flow", Valid = ["on", "off"], Default = "on",Hidden=true);

SET
	Mw = PP.MolecularWeight();
	
	ColumnArea 	= 0.25*Pi*ColumnInternalDiameter^2;
	HETP 		= HeightOfPacking/Number_Stages;	
	V 			= HETP * ColumnArea;
	
	low_pressure = 1E-4 * 'atm';
	low_flow = 1E-6 * 'kmol/h';
	zero_flow = 0 * 'kmol/h';

VARIABLES

	Inlet 			as stream 			(Brief="Feed stream", Symbol="_{in}",Protected=true);
in	InletLiquid 	as stream 			(Brief="Inlet liquid stream",  Symbol="_{inL}",Protected=true);
in	InletVapour 	as stream 			(Brief="Inlet vapour stream",  Symbol="_{inV}",Protected=true);
out	OutletLiquid 	as liquid_stream 	(Brief="Outlet liquid stream", Symbol="_{outL}",Protected=true);
out	OutletVapour 	as vapour_stream 	(Brief="Outlet vapour stream", Symbol="_{outV}",Protected=true);

	M(NComp) 	as mol 			(Brief="Molar Holdup in the tray", Default=0.01, Lower=-0.000001, Upper=100,Protected=true);
	ML 			as mol 			(Brief="Molar liquid holdup", Default=0.01, Lower=0, Upper=100,Protected=true);
	MV 			as mol 			(Brief="Molar vapour holdup", Default=0.01, Lower=0, Upper=100,Protected=true);
	E 			as energy 		(Brief="Total Energy Holdup on tray", Default=-500,Protected=true);
	vL 			as volume_mol 	(Brief="Liquid Molar Volume",Protected=true);
	vV 			as volume_mol 	(Brief="Vapour Molar volume",Protected=true);
	miL 		as viscosity 	(Brief="Liquid dynamic viscosity", DisplayUnit='kg/m/s',Protected=true);
	rhoL 		as dens_mass	(Brief="Liquid mass density",Protected=true);
	rhoV 		as dens_mass	(Brief="Vapour mass density",Protected=true);
	uL 			as velocity 	(Brief="volume flow rate of liquid, m^3/m^2/s", Lower=-10, Upper=1000,Protected=true);
	uV 			as velocity 	(Brief="volume flow rate of vapor, m^3/m^2/s", Lower=-10, Upper=1000,Protected=true);
	Al 			as area 		(Brief="Area occupied by the liquid", Default=0.001, Upper=10,Protected=true);
	hl 			as positive 	(Brief="Column holdup", Unit='m^3/m^3', Default=0.01,Upper=10,Protected=true);
	deltaP 		as pressure 	(Brief="Stage Pressure drop",Protected=true);

EQUATIONS

switch VapourFlow
	
	case "on":
"Pressure drop and Vapor flow, Billet (4-58)"
	deltaP/HETP  = LiquidResistanceCoeff *( 0.5*AreaPerPackingVolume + 2/ColumnInternalDiameter) * 1/((PackingVoidFraction-hl)^3) * (uV^2) *rhoV;
	
	when InletVapour.F < low_flow switchto "off";
	
	case "off":
"Vapour Flow"
	InletVapour.F = zero_flow;
	
	when deltaP > low_pressure switchto "on";

end

"Energy Balance"
	diff(E) = (Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h- OutletLiquid.F*OutletLiquid.h
	- OutletVapour.F*OutletVapour.h + HeatOnStage );

"Energy Holdup"
	E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletLiquid.P*V;

"Geometry Constraint"
	V*PackingVoidFraction= ML*vL + MV*vV;

"Volume flow rate of vapor, m^3/m^2/s"
	uV * (V*PackingVoidFraction/HETP - Al) = InletVapour.F * vV;

"Liquid holdup"
	hl*V*PackingVoidFraction = ML*vL;

"Liquid velocity as a function of liquid holdup, Billet (4-27)"
	hl^3 = (12/Gconst) * AreaPerPackingVolume^2 * (miL/rhoL) * uL;

"Area occupied by the liquid"
	Al = ML*vL/HETP;

"Component Molar Balance"
	diff(M)=Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z;

"Molar Holdup"
	M = ML*OutletLiquid.z + MV*OutletVapour.z;

"Mol Fraction Normalisation"
	sum(OutletLiquid.z)= 1.0;

"Mol Fraction Constraint"
	sum(OutletLiquid.z)= sum(OutletVapour.z);
	
"Liquid Volume"
	vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Vapour Volume"
	vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z);
	
"Chemical Equilibrium"
	PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = 	PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, OutletVapour.z)*OutletVapour.z;
	
"Thermal Equilibrium"
	OutletVapour.T = OutletLiquid.T;
	
"Mechanical Equilibrium"
	OutletLiquid.P = OutletVapour.P;
	
"Liquid Density"
	rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Vapour Density"
	rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z);

"Liquid viscosity"
	miL = PP.LiquidViscosity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Volume flow rate of liquid, m^3/m^2/s"
	uL * Al = OutletLiquid.F * vL;
	
"Pressure Drop"
	deltaP = InletVapour.P - OutletVapour.P;	

end

#*-------------------------------------
* Nonequilibrium Model
-------------------------------------*
Model interfaceTeste
	
	ATTRIBUTES
	Pallete 	= false;
	Icon 		= "icon/Tray"; 
	Brief 		= "Descrition of variables of the equilibrium interface.";
	Info 		=
"This model contains only the variables of the equilibrium interface.";

	PARAMETERS
outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
outer NComp as Integer;
outer NC1 as Integer;
	
	VARIABLES
	NL(NComp) as flow_mol_delta	(Brief = "Stream Molar Rate on Liquid Phase");
	NV(NComp) as flow_mol_delta	(Brief = "Stream Molar Rate on Vapour Phase");
	T as temperature 		(Brief = "Stream Temperature");
	P as pressure 			(Brief = "Stream Pressure");
	x(NComp) as fraction	(Brief = "Stream Molar Fraction on Liquid Phase");
	y(NComp) as fraction	(Brief = "Stream Molar Fraction on Vapour Phase");
	a as area 			    (Brief = "Interface Area");
	htL as heat_trans_coeff (Brief = "Heat Transference Coefficient on Liquid Phase");
	htV as heat_trans_coeff (Brief = "Heat Transference Coefficient on Vapour Phase");	
	E_liq as heat_rate      (Brief = "Liquid Energy Rate at interface");
    E_vap as heat_rate      (Brief = "Vapour Energy Rate at interface");	
	hL as enth_mol 	       	(Brief = "Liquid Molar Enthalpy");
	hV as enth_mol 	       	(Brief = "Vapour Molar Enthalpy");
	kL(NC1,NC1) as velocity (Brief = "Mass Transfer Coefficients");
	kV(NC1,NC1) as velocity (Brief = "Mass Transfer Coefficients");
	
	EQUATIONS
	"Liquid Enthalpy"
	hL = PP.LiquidEnthalpy(T, P, x);
	
	"Vapour Enthalpy"
	hV = PP.VapourEnthalpy(T, P, y);

end

Model trayRateBasicTeste
	ATTRIBUTES
	Pallete 	= false;
	Icon 		= "icon/Tray"; 
	Brief 		= "Basic equations of a tray rate column model.";
	Info 		=
"This model contains only the main equations of a column tray nonequilibrium model without
the hidraulic equations.
	
== Assumptions ==
* both phases (liquid and vapour) exists all the time;
* no entrainment of liquid or vapour phase;
* no weeping;
* the dymanics in the downcomer are neglected.
";
	
	PARAMETERS
outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
outer NComp as Integer;
    NC1 as Integer;
	V as volume(Brief="Total Volume of the tray");
	Q as heat_rate (Brief="Rate of heat supply"); 
	Ap as area (Brief="Plate area = Atray - Adowncomer");
	
	VARIABLES
in	Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
in	InletFV as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
in	InletLiquid as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
in	InletVapour as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
out	OutletLiquid as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
out	OutletVapour as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");

	M_liq(NComp) as mol (Brief="Liquid Molar Holdup in the tray");
	M_vap(NComp) as mol (Brief="Vapour Molar Holdup in the tray");
	ML as mol (Brief="Molar liquid holdup");
	MV as mol (Brief="Molar vapour holdup");
	E_liq as energy (Brief="Total Liquid Energy Holdup on tray");
	E_vap as energy (Brief="Total Vapour Energy Holdup on tray");
	vL as volume_mol (Brief="Liquid Molar Volume");
	vV as volume_mol (Brief="Vapour Molar volume");
	Level as length (Brief="Height of clear liquid on plate");
	interf as interfaceTeste;	

	SET   
	NC1=NComp-1;

	EQUATIONS
	"Component Molar Balance"
	diff(M_liq)=Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z 
	- OutletLiquid.F*OutletLiquid.z + interf.NL;
	
	diff(M_vap)=InletFV.F*InletFV.z + InletVapour.F*InletVapour.z 
	- OutletVapour.F*OutletVapour.z - interf.NV;
	
	"Energy Balance"
	diff(E_liq) = Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h
		- OutletLiquid.F*OutletLiquid.h  + Q + interf.E_liq;
	
	diff(E_vap) = InletFV.F*InletFV.h + InletVapour.F*InletVapour.h
		- OutletVapour.F*OutletVapour.h  - interf.E_vap;
	
	"Molar Holdup"
	M_liq = ML*OutletLiquid.z;
	
	M_vap = MV*OutletVapour.z;
	
	"Energy Holdup"
	E_liq = ML*(OutletLiquid.h - OutletLiquid.P*vL);
	
	E_vap = MV*(OutletVapour.h - OutletVapour.P*vV);
	
	"Energy Rate through the interface"
	interf.E_liq = interf.htL*interf.a*(interf.T-OutletLiquid.T)+sum(interf.NL)*interf.hL;	
	
	interf.E_vap = interf.htV*interf.a*(OutletVapour.T-interf.T)+sum(interf.NV)*interf.hV;
	
	"Mass Conservation"
	interf.NL = interf.NV;
	
	"Energy Conservation"
	interf.E_liq = interf.E_vap;
	
	"Mol fraction normalisation"
	sum(OutletLiquid.z)= 1.0;
	sum(OutletLiquid.z)= sum(OutletVapour.z);
	sum(interf.x)=1.0;
	sum(interf.x)=sum(interf.y);
	
	"Liquid Volume"
	vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);
	"Vapour Volume"
	vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z);
	
	"Chemical Equilibrium"
	PP.LiquidFugacityCoefficient(interf.T, interf.P, interf.x)*interf.x = 
		PP.VapourFugacityCoefficient(interf.T, interf.P, interf.y)*interf.y;

	"Geometry Constraint"
	V = ML*vL + MV*vV;
	
	"Level of clear liquid over the weir"
	Level = ML*vL/Ap;

	"Total Mass Transfer Rates"
	interf.NL(1:NC1)=interf.a*sumt(interf.kL*(interf.x(1:NC1)-OutletLiquid.z(1:NC1)))/vL+
		OutletLiquid.z(1:NC1)*sum(interf.NL);

#	interf.NL(1:NC1)=0.01*'kmol/s';
	
	interf.NV(1:NC1)=interf.a*sumt(interf.kV*(OutletVapour.z(1:NC1)-interf.y(1:NC1)))/vV+
		OutletVapour.z(1:NC1)*sum(interf.NV);

	"Mechanical Equilibrium"
	OutletVapour.P = OutletLiquid.P;
	interf.P=OutletLiquid.P;
end

Model trayRateTeste as trayRateBasicTeste
	ATTRIBUTES
	Pallete 	= false;
	Icon 		= "icon/Tray"; 
	Brief 		= "Complete rate model of a column tray.";
	Info 		=
"== Specify ==
* the Feed stream
* the Liquid inlet stream
* the Vapour inlet stream
* the Vapour outlet flow (OutletVapour.F)
	
== Initial ==
* the plate temperature of both phases (OutletLiquid.T and OutletVapour.T)
* the liquid height (Level) OR the liquid flow holdup (ML)
* the vapor holdup (MV)
* (NoComps - 1) OutletLiquid compositions
";

	PARAMETERS
	Ah as area (Brief="Total holes area");
	lw as length (Brief="Weir length");
	g as acceleration (Default=9.81);
	hw as length (Brief="Weir height");
	beta as fraction (Brief="Aeration fraction");
	alfa as fraction (Brief="Dry pressure drop coefficient");
	
	VapourFlow as Switcher(Valid = ["on", "off"], Default = "on");
	LiquidFlow as Switcher(Valid = ["on", "off"], Default = "on");
	
	VARIABLES
	rhoL as dens_mass;
	rhoV as dens_mass;

	EQUATIONS
	"Liquid Density"
	rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);
	"Vapour Density"
	rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z);

	switch LiquidFlow
		case "on":
		"Francis Equation"
#		OutletLiquid.F*vL = 1.84*'m^0.5/s'*lw*((Level-(beta*hw))/(beta))^1.5;
		OutletLiquid.F*vL = 1.84*'1/s'*lw*((Level-(beta*hw))/(beta))^2;
		when Level < (beta * hw) switchto "off";
		
		case "off":
		"Low level"
		OutletLiquid.F = 0 * 'mol/h';
		when Level > (beta * hw) + 1e-6*'m' switchto "on";
	end

	switch VapourFlow
		case "on":
		InletVapour.F*vV = sqrt((InletVapour.P - OutletVapour.P)/(rhoV*alfa))*Ah;
		when InletVapour.F < 1e-6 * 'kmol/h' switchto "off";
		
		case "off":
		InletVapour.F = 0 * 'mol/s';
		when InletVapour.P > OutletVapour.P + Level*g*rhoL + 1e-1 * 'atm' switchto "on";
	end	
end

*#