#*---------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* 9. Gas phase catalytic reactor
*----------------------------------------------------------------------
*
*   Description:
*      This problem is part of a collection of 10 representative
*	problems in Chemical Engineering for solution by numerical methods
*	developed for Cutlip (1998).
*
*   Subject:
*       * Reaction Engineering
*
*	Concepts utilized:
*		Design of a gas phase catalytic reactor with pressure drop for
*	a first order reversible gas phase reaction.
*
*	Numerical method:
*		* Simultaneous ODEs with known boundary conditions
*
*   Reference:
*		* CUTLIP et al. A collection of 10 numerical problems in 
*	chemical engineering solved by various mathematical software 
*	packages. Comp. Appl. in Eng. Education. v. 6, 169-180, 1998. 
*       * More informations and a detailed description of all problems
*	is available online in http://www.polymath-software.com/ASEE
* 
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* GIMSCOP/UFRGS - Group of Integration, Modeling, Simulation, 
*					Control, and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#
using "types";



Model stream
	PARAMETERS
outer NComp 	as Integer		(Brief="Number of components", Lower=1);
	cp(NComp)	as cp_mol;
	
	SET
	cp = [40, 80]*'J/mol/K';
	
	VARIABLES
	C(NComp)	as conc_mol		(Brief="Molar concentration", Lower=0, DisplayUnit='mol/l');
	T			as temperature;
	P			as pressure;
end



FlowSheet reactor
	PARAMETERS
	NComp 		as Integer		(Brief="Number of components", Lower=1);
	NReac		as Integer		(Brief="Number of reactions");
	
	Tref		as temperature	(Brief="Reference temperature", Default=450);
	Ta			as temperature;
	
	ko(NReac)	as Real 		(Brief="Frequency factor at Treff", Unit='l^2/kg/min/mol');
	Ko(NReac)	as Real			(Unit='l/mol');
	E(NReac)	as energy_mol 	(Brief="Activation energy");
	R			as Real			(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
	
	DH(NReac)	as enth_mol;
	U			as Real			(Unit='J/kg/K/min');
	alpha		as Real			(Unit='1/kg');
	
	
	VARIABLES
	Inlet		as stream;
	Outlet		as stream;
	
	Fo(NComp)	as flow_mol		(Brief="Inlet molar flow");
	X			as fraction		(Brief="Molar conversion of A");
	y			as fraction		(Brief="Normalized pressure");
	Tn			as Real			(Brief="Temperature by 1000", Lower=0.273);
	
	
	r(NReac)	as Real			(Brief="Mass rate of reaction", Unit='mol/kg/min');
	k(NReac)	as Real			(Brief="Specific rate of reaction", Unit='l^2/kg/min/mol');
	K(NReac)	as Real			(Brief="Equilibrium constant", Unit='l/mol');
	
	W			as mass			(Brief="Catalytic weight");
	
	
	EQUATIONS
	"Change time in W"
	W = time*'kg/s';
	
	"Mole balance"
	diff(X) = -r(1)/Fo(1)*'kg/s';
	
	"Rate of reaction"
	-r = k*(Outlet.C(1)^2 - Outlet.C(2)/K);
	
	
	"Specific rate of reaction"
	k = ko*exp(E/R*(1/Tref - 1/Outlet.T));
	
	"Equilibrium constant"
	K = Ko*exp(DH/R*(1/Tref - 1/Outlet.T));
	
	
	"Molar concentration of A"
	Outlet.C(1) = Inlet.C(1)*(1 - X)/(1 - 0.5*X)*y*(Inlet.T/Outlet.T);
	
	"Molar concentration of C"
	Outlet.C(2) = 0.5*Inlet.C(1)*X/(1 - 0.5*X)*y*(Inlet.T/Outlet.T);
	
	
	"Pressure drop"
	diff(y) = -alpha*(1 - 0.5*X)/2/y*(Outlet.T/Inlet.T)*'kg/s';
	
	"Energy balance"
	diff(Outlet.T) = (U*(Ta - Outlet.T) + r*DH)/(Fo(1)*Inlet.cp(1))*'kg/s';
	
	
	"Normalized pressure"
	y = Outlet.P/Inlet.P;
	
	"Temperature/1000"
	Tn = Outlet.T/1e3/'K';
	
	
	SET
	NComp = 2; # A, and C
	
	DH(1) = -4e4*'J/mol';
	ko(1) = 0.5*'l^2/kg/min/mol'; # at Tref
	Ko(1) = 2.5e4*'l/mol'; # at Tref
	E(1)  = 4.18e4*'J/mol';
	
	Ta = 500*'K';
	U = 0.8*'J/kg/K/min';
	alpha = 0.015/'kg';

	
	SPECIFY
	Fo = [5, 0]*'mol/min';
	Inlet.C = [0.271, 0]*'mol/l'; 
	Inlet.T = 450*'K';
	Inlet.P = 10*'atm';

	
	INITIAL
	"Molar conversion"
	X = 0;
	
	"Drop pressure"
	y = 1;
	
	"Temperature"
	Outlet.T = Inlet.T;


	OPTIONS
	TimeStart = 0;
	TimeStep = 0.25;
	TimeEnd = 20;
	
	TimeUnit = 's';
end