#*-------------------------------------------------------------------
* 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.
*
*---------------------------------------------------------------------
* Production of acetic anhydride
*----------------------------------------------------------------------
* Solved problem from Fogler (1999)
* Problem number: 8-7
* Page: 421 (Brazilian version, 2002)
*----------------------------------------------------------------------
*
*   Description:
*		The acetic anhydride is produced for thermal craking of the
*	acetone in a PFR:
*			CH3COCH3 -> CH2CO + CH4
*		This sample calculates the molar conversion and temperature
*	as function of the length in the tubular reactor. In the case I 
*	the operation is adiabatic and in the case II the reactor is
*	jacketed.
*
*   Assumptions
*		* first-order reaction with respect to acetone
*		* steady-state
*       * gaseous phase
*
*	Specify:
*		* the inlet stream
*		* the kinetic parameters
*		* the parameters of components
*
*----------------------------------------------------------------------
* Author: Christiano D. W. Guerra and Rodolfo Rodrigues
* $Id: acetic_anhydride.mso 574 2008-07-25 14:18:50Z rafael $
*--------------------------------------------------------------------*#

using "types";


#*---------------------------------------------------------------------
* Model of the thermal craking of acetone
*--------------------------------------------------------------------*#

Model thermal_cracking
	PARAMETERS
	NComp		as Integer		(Brief="Number of components", Lower=1);
	stoic(NComp)as Real			(Brief="Stoichiometric number");
	Pao 		as pressure 	(Brief="Input pressure of A");
	Tr			as temperature	(Brief="Reference temperature");
	To			as temperature	(Brief="Inlet temperature");
	Ta			as temperature	(Brief="Internal temperature");
	Hr(NComp)	as enth_mol		(Brief="Enthalpy of component");
	R			as Real			(Brief="Universal gas constant", Unit='kPa*m^3/kmol/K', Default=8.314);
	U			as Real			(Brief="Heat transfer coefficient", Unit='J/m^2/K/s');
	a			as Real			(Brief="Heat transfer area per volume of tube", Unit='1/m');
	alpha(NComp)as cp_mol		(Brief="Alpha term of Cp expression");
	beta(NComp)	as Real			(Brief="Beta term of Cp expression", Unit='J/mol/K^2');
	gamma(NComp)as Real			(Brief="Gamma term of Cp expression", Unit='J/mol/K^3');
	
	VARIABLES
	Ca			as conc_mol 	(Brief="Molar concentration of A", DisplayUnit='kmol/m^3');
	Cao 		as conc_mol 	(Brief="Inlet molar concentration of A", DisplayUnit='mol/m^3');
	Fao 		as flow_mol 	(Brief="Inlet molar flow of A");
	vo			as flow_vol   	(Brief="Volumetric flow", DisplayUnit='m^3/s');
	r 			as reaction_mol	(Brief="Rate of reaction", DisplayUnit='kmol/m^3/s');
	k 			as Real			(Brief="Specific rate of reaction", Unit='1/s');
	T 			as temperature 	(Brief="Temperature of reactor");
	X 			as fraction 	(Brief="Molar conversion");
	V			as volume		(Brief="Volume");
	eps			as Real			(Brief="Parameter epsilon");
	Cp(NComp)	as cp_mol		(Brief="Molar heat capacity", DisplayUnit='J/mol/K');
	DHr			as enth_mol		(Brief="Enthalpy of reaction", DisplayUnit='kJ/mol');	
	
	EQUATIONS
	"Change time in V"
	V = time*'m^3/s';
	
	"Molar balance"
	diff(X) = (-r)/Fao*'m^3/s';
	
	"Rate of reaction"
	r = -k*Ca;
	
	"Specific rate of reaction"
	k = exp(34.34)*exp(-34222*'K'/T)*'1/s';
	
	"Concentration of component A"
	Ca = Cao*(1 - X)/(1 + eps*X)*To/T;
	
	"Parameter epsilon"
	eps = sum(stoic); # yAo = 1
	
	"Inlet molar concentration of A"
	Cao = Pao/(R*To);
	
	"Volumetric flow"
	Fao = Cao*vo;
	
	"Energy balance"
	diff(T)*(Fao*(Cp(1) + X*sumt(stoic*Cp))) = (U*a*(Ta - T) + (-r)*(-DHr))*'m^3/s';

	"Enthalpy of reaction"
	DHr = sumt(stoic*Hr) + sumt(stoic*alpha)*(T - Tr) + sumt(stoic*beta)/2*(T^2 - Tr^2) 
		+  sumt(stoic*gamma)/3*(T^3 - Tr^3);
	
	"Molar heat capacity"
	Cp  = alpha + beta*T + gamma*T^2;
end


#*---------------------------------------------------------------------
* Case I: In an adiabatic PFR
*--------------------------------------------------------------------*#

FlowSheet adiabatic_reactor
	DEVICES
	R as thermal_cracking;
	
	SET
	R.NComp = 3; # A: acetone, B: ketene and C: methane
	R.stoic = [-1.0, 1.0, 1.0]; # A -> B + C
	
	R.Pao = 162*'kPa';
	
	R.alpha = [26.63,  20.04, 13.39]*'J/mol/K';
	R.beta  = [0.183, 0.0945, 0.077]*'J/mol/K^2';
	R.gamma = [-45.86e-6, -30.95e-6, -18.71e-6]*'J/mol/K^3';
	
	R.Hr = [-216.67, -61.09, -74.81]*'kJ/mol';
	R.Tr = 298*'K';
	R.To = 1035*'K';
	R.Ta = 1150*'K';
	R.U  = 0.0*'J/m^2/K/s';
	R.a  = 150*'1/m';
	
	SPECIFY
	R.Fao = (8000/58)*'kmol/h';
	
	INITIAL
	"Molar conversion"
	R.X = 0.0;
	"Temperature"
	R.T = 1035*'K';
	
	OPTIONS
	TimeStep = 0.05;
	TimeEnd = 5;
	TimeUnit = 's';
end


#*---------------------------------------------------------------------
* Case II: In an jacketed PFR
*--------------------------------------------------------------------*#

FlowSheet jacketed_reactor
	DEVICES
	R as thermal_cracking;
	
	SET
	R.NComp = 3; # A: acetone, B: ketene and C: methane
	R.stoic = [-1.0, 1.0, 1.0]; # A -> B + C
	
	R.Pao = 162*'kPa';
	
	R.alpha = [26.63,  20.04, 13.39]*'J/mol/K';
	R.beta  = [0.183, 0.0945, 0.077]*'J/mol/K^2';
	R.gamma = [-45.86e-6, -30.95e-6, -18.71e-6]*'J/mol/K^3';
	
	R.Hr = [-216.67, -61.09, -74.81]*'kJ/mol';
	R.Tr = 298*'K';
	R.To = 1035*'K';
	R.Ta = 1150*'K';
	R.U  = 110*'J/m^2/K/s';
	R.a  = 150*'1/m';
	
	SPECIFY
	R.Fao = (18.8*2e-3)*'mol/s';
	
	INITIAL
	"Molar conversion"
	R.X = 0.0;
	"Temperature"
	R.T = 1035*'K';
	
	OPTIONS
	TimeStep = 1e-5;
	TimeEnd = 1e-3;
	TimeUnit = 's';
end