#*---------------------------------------------------------------------
* 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 an equilibrium reactor
*----------------------------------------------------------------------
*
*   Description:
*       Thermodynamic equilibrium modeling of a reactor based on 
*	equilibrium constants approach.
*
*   Assumptions:
*       * thermodynamic equilibrium
*		* steady-state
*
*	Specify:
*		* inlet stream
*		* stoichiometric matrix
*		* equilibrium temperature
*
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* $Id$
*--------------------------------------------------------------------*#

using "tank_basic";


#*---------------------------------------------------------------------
*	only vapour phase
*--------------------------------------------------------------------*#
Model equil_vap as tank_vap
	ATTRIBUTES
	Pallete 	= true;
	Icon 		= "icon/cstr"; 
	Brief 		= "Model of a generic vapour-phase equilibrium CSTR";
	Info 		= "
Requires the information of:
* number of reactions
* matrix of stoichiometric coefficients (components by reactions)
";

	PARAMETERS
	NReac		as Integer 		(Brief="Number of reactions", Default=1);
    stoic(NComp,NReac) as Real 	(Brief="Stoichiometric matrix");
	Rg			as Real			(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
	fs(NComp) 	as pressure 	(Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
	To			as temperature  (Brief="Reference temperature", Default=298.15);

	VARIABLES
out Outlet		as vapour_stream; # Outlet stream

	G(NComp)	as energy_mol 	(Brief="Gibbs free-energy of formation");
	K(NReac)	as Real			(Brief="Equillibrium constant",Default=1.5);
	activ(NComp)as Real 		(Brief="Activity",Default=0.2);
	
	rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
	extent(NReac) as flow_mol 	(Brief="Extent of reaction");
	conv(NComp) as Real 		(Brief="Fractional conversion of component", Default=0); # Lower=-1e3, Upper=1e3);
	
	EQUATIONS
	"Outlet stream"
	Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
	
	"Mechanical equilibrium"
	Outlet.P = Outletm.P;
	
	"Energy balance"
	Outlet.F*Outlet.h = Outletm.F*Outletm.h;
	
	"Steady-state"
	Outlet.F = Inlet.F + sum(sumt(stoic*extent));

	"Gibbs free-energy of formation"
	G = PP.IdealGasGibbsOfFormation(Outlet.T);

#	"Gibbs free-energy of formation without Cp correction"
#	G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);

	"Gibbs free energy of reaction"
#	sumt(G*stoic) = -Rg*Outlet.T*ln(K);
	K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
	
	for j in [1:NReac]
	"Equilibrium constant"
		K(j) = prod(activ^stoic(:,j));
	end

	for i in [1:NComp]
	"Outlet molar fraction"
		Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
	end	

	for i in [1:NComp]
	  if (Outletm.z(i) > 0) then
	    "Molar conversion"
	    Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
	  else if (Outlet.z(i) > 0) then
			"Molar conversion"
				conv(i) = 1;	# ?
			else
			"Molar conversion"
				conv(i) = 0;	# ?
			end
	  end
	end

	"Activity"
	activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
end


#*---------------------------------------------------------------------
*	only liquid-phase
*--------------------------------------------------------------------*#
Model equil_liq as tank_liq
	ATTRIBUTES
	Pallete 	= true;
	Icon 		= "icon/cstr"; 
	Brief 		= "Model of a generic liquid-phase equilibrium CSTR";
	Info 		= "
Requires the information of:
* number of reactions
* matrix of stoichiometric coefficients (components by reactions)
";

	PARAMETERS
	NReac		as Integer 		(Brief="Number of reactions", Default=1);
    stoic(NComp,NReac) as Real 	(Brief="Stoichiometric matrix");
	Rg			as Real			(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
	Ps			as pressure 	(Brief="Standard pressure", Default=1, DisplayUnit='bar');
	To			as temperature  (Brief="Reference temperature", Default=298.15);

	VARIABLES
out Outlet		as liquid_stream; # Outlet stream

	G(NReac)	as enth_mol 	(Brief="Gibbs free-energy of formation");
	K(NReac)	as fraction		(Brief="Equillibrium constant");
	activ(NComp)as Real 		(Brief="Activity");
	
	rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
	extent(NReac)as flow_mol 	(Brief="Extent of reaction");
	conv(NComp) as Real 		(Brief="Fractional conversion of component", Default=0);
	
	EQUATIONS
	"Outlet stream"
	Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
	
	"Mechanical equilibrium"
	Outlet.P = Outletm.P;
	
	"Energy balance"
	Outlet.F*Outlet.h = Outletm.F*Outletm.h;
	
	"Steady-state"
	Outlet.F = Inlet.F + sum(sumt(stoic*extent));
	
	"Gibbs free-energy of formation"
	G = PP.IdealGasGibbsOfFormation(Outlet.T);

#	"Gibbs free-energy of formation without Cp correction"
#	G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);

	"Gibbs free energy of reaction"
#	sumt(G*stoic) = -Rg*Outlet.T*ln(K);
	K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
	
	for j in [1:NReac]
	"Equilibrium constant"
		K(j) = prod(activ^stoic(:,j));
	end
	
	for i in [1:NComp]
	"Outlet molar fraction"
		Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
	end	

	for i in [1:NComp]
	  if (Outletm.z(i) > 0) then
	    "Molar conversion"
	    Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
	  else if (Outlet.z(i) > 0) then
			"Molar conversion"
				conv(i) = 1;	# ?
			else
			"Molar conversion"
				conv(i) = 0;	# ?
			end
	  end
	end
	
	"Activity"
	activ = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z
		*exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
end