#*---------------------------------------------------------------------
* 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 simplified Gibbs reactor.
* This model requires VRTherm (www.vrtech.com.br) to run.
*----------------------------------------------------------------------
*
*
*
*----------------------------------------------------------------------
* Author: Rafael de Pelegrini Soares
* $Id$
*--------------------------------------------------------------------*#

using "types";

Model gibbs_reactor_simple
	ATTRIBUTES
	Info =
"Model for reactor in thermodynamic equilibrium based on Gibbs free
energy.

The user should specify T, P and ni.

There is a correction for G0 as a function of the temperature.
This correction considers that H0 does not depend on the temperature.
This is a good approximation for most cases (less than 2% of error)";
	
	PARAMETERS
outer PP		as Plugin (Brief="External physical properties", Type="PP");
outer NComp 	as Integer (Brief="Number of components", Default=1);

	nu(NComp) as Real(Symbol="\nu_i");

	R  as Real(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
	T0 as temperature(Default = 298.15, Symbol="T_{298}");
	P0 as Real(Default=1, Unit='bar', Symbol="P^0");
	g0(NComp) as energy_mol 	(Brief="Gibbs energy in standard state", Symbol="g^0_{298,i}");
	h0(NComp) as energy_mol 	(Brief="Enthalpy in standard state", Symbol="h^0_{298,i}");
	
	SET
	g0 = PP.IdealGasGibbsOfFormationAt25C();
	h0 = PP.IdealGasEnthalpyOfFormationAt25C();

	VARIABLES
	T as temperature;
	P as pressure;
	n0(NComp) as positive(Brief="Initial number of mols", Unit='mol', Symbol="n_{i,0}");
	n(NComp)  as positive(Brief="Number of mols at equilibrium", Unit='mol', Symbol="n_i");
	advance   as Real(Unit='mol', Symbol="\epsilon");

	K  as positive(Brief="Reaction equilibrium constant at T");
	K0 as positive(Brief="Reaction equilibrium constant at T0", Symbol="K_{298}");
	
	phi(NComp)	as fugacity(Brief="Fugacity coefficient", Default=1, Symbol="\phi_i");
	
	EQUATIONS
	"Equilibrium constant at 298 K"
	K0 = exp(-sum(nu*g0)/(R*T0));
	"Equilibrium constant at T"
	K = K0 * exp(-sum(nu*h0)/R*(1/T - 1/T0));

	"Equilibrium rule"
	K = prod( (n/sum(n)*phi*P/P0) ^ nu);

	"Reaction advance"
	n = n0 + nu * advance;
	
	"Fugacity coefficient"
	phi = PP.VapourFugacityCoefficient(T, P, n/sum(n));
end

# Ethane decomposition to produce ethylene and hydrogen at 1000 degC
FlowSheet gibbs_reactor_simple_sample
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["ethane", "ethylene", "hydrogen"],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;

	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, 1, 1];

	SPECIFY
	R.T = (1000 + 273.15) * 'K';
	R.P = 1 * 'atm';
	R.n0 = [1, 0, 0] * 'mol';
	
	# Expected results:
	# advance = 0.9;
	# n = [0.1, 0.9, 0.9]
	
	OPTIONS
	Dynamic = false;
end


# Ethanol production by ethylene hydratation at 250 degC and 35 bar
FlowSheet gibbs_reactor_simple_sample2
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["ethylene", "water", "ethanol"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;

	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -1, 1];

	SPECIFY
	R.T = (250 + 273.15) * 'K';
	R.P = 35 * 'atm';
	R.n0 = [1, 5, 0] * 'mol';
	
	# Expected results:
	# advance = 0.21;
	
	OPTIONS
	Dynamic = false;
end

# Water-gas-shift T = 1100 K, P = 1 bar
FlowSheet gibbs_reactor_simple_sample3
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;

	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -1, 1, 1];

	SPECIFY
	R.T = 1100 * 'K';
	R.P = 1 * 'bar';
	R.n0 = [1, 1, 0, 0] * 'mol';
	
	# Expected results:
	# advance = 0.41;
	
	OPTIONS
	Dynamic = false;
end


# Water-gas-shift T = 1100 K, P = 10 bar
FlowSheet gibbs_reactor_simple_sample4
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;

	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -1, 1, 1];

	SPECIFY
	R.T = 1100 * 'K';
	R.P = 1 * 'bar';
	R.n0 = [1, 1, 0, 0] * 'mol';
	
	# Expected results:
	# advance = 0.41;
	
	OPTIONS
	Dynamic = false;
end


# Water-gas-shift T = 1100 K, P = 1 bar excess of water
FlowSheet gibbs_reactor_simple_sample5
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;
	
	DEVICES
	R as gibbs_reactor_simple;
	
	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -1, 1, 1];

	SPECIFY
	R.T = 1100 * 'K';
	R.P = 1 * 'bar';
	R.n0 = [1, 2, 0, 0] * 'mol';
	
	# Expected results:
	# advance = 0.56;
	
	OPTIONS
	Dynamic = false;
end

# Ammonia synthesis from nitrogen and hydrogen T = 500 degC, P = 1 bar
FlowSheet gibbs_reactor_simple_sample6
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["nitrogen", "hydrogen", "ammonia"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;
	
	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -1, 2];

	SPECIFY
	R.T = (500 + 273.15) * 'K';
	R.P = 1 * 'bar';
	R.n0 = [1, 3, 0] * 'mol';
	
	# Expected results:
	# K = 6e-5;
	
	OPTIONS
	Dynamic = false;
end

# Ammonia synthesis from nitrogen and hydrogen T = 500 degC, P = 300 bar
FlowSheet gibbs_reactor_simple_sample7
	PARAMETERS
	PP as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["nitrogen", "hydrogen", "ammonia"],
		LiquidModel = "IdealLiquid",
		VapourModel = "PR"
	);
	NComp as Integer;

	DEVICES
	R as gibbs_reactor_simple;
	
	SET
	NComp = PP.NumberOfComponents;
	R.nu = [-1, -3, 2];

	SPECIFY
	R.T = (500 + 273.15) * 'K';
	R.P = 300 * 'bar';
	R.n0 = [1, 3, 0] * 'mol';
	
	# Expected results:
	# advance = 0.54;
	
	OPTIONS
	Dynamic = false;
end