#*---------------------------------------------------------------------
* 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 Gibbs reactor
*----------------------------------------------------------------------
*
*
*
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* $Id$
*--------------------------------------------------------------------*#

using "tank_basic";


#*---------------------------------------------------------------------
*	only vapour phase
*--------------------------------------------------------------------*#

Model gibbs_vap as tank_vap
	ATTRIBUTES
	Pallete 	= true;
	Icon 		= "icon/cstr"; 
	Brief 		= "Model of a generic vapour-phase Gibbs CSTR";
	Info 		= "
Requires the information of:
* number of elements
* matrix of elements (elements, compounds)
";

	PARAMETERS
outer NElem			as Integer 		(Brief="Number of elements", Default=1);
	Rg				as Real			(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
	na(NElem,NComp)	as Real			(Brief="Number of elements per component");
	fs(NComp) 		as pressure 	(Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
	T0				as temperature  (Brief="Reference temperature", Default=298.15);
	
	VARIABLES
out Outlet			as vapour_stream; # Outlet stream

	Gs(NComp)		as energy_mol 	(Brief="Gibbs energy in standard state");
	lambda(NElem)	as energy_mol 	(Brief="Lagrangian multiplier");
	phi(NComp)		as fugacity		(Brief="Fugacity coefficient", Default=1);
	activ(NComp)	as Real 		(Brief="Activity", Lower=1e-20);

	rate(NComp) 	as reaction_mol (Brief="Overall component rate of reaction");
	conv(NComp)   	as Real 		(Brief="Fractional conversion of component", Default=0);
	Fi(NComp)		as flow_mol		(Brief="Component molar flow rate");

	EQUATIONS
	"Outlet stream"
	Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
	
	"Mechanical equilibrium"
	Outlet.P = Outletm.P;

	"Steady-state"
	Outlet.F = sum(Fi);

	"Component molar flow rate"
	Fi = Outlet.F*Outlet.z;
	
	"Energy balance"
	Outlet.F*Outlet.h = Outletm.F*Outletm.h;

	"Element balance"
	sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);

	"Gibbs Energy of Formation"
	Gs = PP.IdealGasGibbsOfFormation(Outlet.T);

#	"Gibbs Energy of Formation without Cp correction"
#	Gs = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/T0+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/T0);

	for i in [1:NComp]
	"Lagrangian multiplier"
		Gs(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));

	  if (Outletm.z(i) > 0) then
	    "Molar conversion"
	    Fi(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
	
	"Fugacity coefficient"
	phi = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
	
	"Activity"
	activ = phi*Outlet.P*Outlet.z/fs;
end