#*-------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* Author: Rafael de P. Soares and Paula B. Staudt
* $Id: pfr.mso 745 2009-03-20 19:59:07Z bicca $
*--------------------------------------------------------------------*#

using "streams";

Model pfr

ATTRIBUTES
	Pallete 	= true;
	Brief 		= "Model of a Generic PFR with constant mass holdup";
	Icon		= "icon/pfr";
	Info 		= 
"== Requires the information of ==
* Reaction values
* Heat of reaction
* Pressure profile
";

PARAMETERS

outer 	PP   							as Plugin 		(Brief = "External Physical Properties", Type="PP");
outer	NComp 						as Integer		(Brief="Number of components");
		
	NReac 							as Integer		(Brief="Number of reactions");
	stoic(NComp, NReac) 	as Real 			(Brief = "Stoichiometric Matrix");
	NDisc 								as Integer		(Brief="Number of points of discretization", Default=10);
	Mw(NComp)  					as molweight 	(Brief="Component Mol Weight");
	L 										as length			(Brief="Reactor Length");
	Across 							as area			(Brief="Cross section area");
	Dpipe 								as length			(Brief="Reactor Inner Diameter");
	pi   									as Real 			(Brief="pi number",Default=3.141592, Symbol = "\pi");
	dx 									as length 		(Brief = "Incremental Length");
	dv 									as volume 		(Brief = "Incremental volume");
	Roughness 					as length			(Brief="Reactor Tube Roughness", Default = 4.572E-5, Symbol = "\varepsilon");

SET
	
	Mw   		= PP.MolecularWeight();
	Across 	= 0.25*pi*Dpipe^2; 
	dx			= L/NDisc;
	dv			= Across*dx;

VARIABLES

in		Inlet  	as stream	(Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
out	Outlet 	as stream	(Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}");

	str(NDisc+1) 			as streamPH;
	vel(NDisc+1) 			as velocity;
	vol(NDisc+1) 			as vol_mol;
	rho(NDisc+1) 			as dens_mass;
	rhom(NDisc+1) 		as dens_mol;
	dP(NDisc+1) 			as press_delta 		(Brief = "Friction Pressure Drop");
	Lincr(NDisc+1) 		as length 			(Brief = "Length Points", Symbol = "L_{incr}");
	
	q(NDisc)    						as heat_rate;
	C(NComp, NDisc+1)  	as conc_mol			(Brief="Components concentration");
	E(NDisc) 						as energy 				(Brief="Total Energy Holdup on element");
	r(NReac, NDisc) 			as reaction_mol;
	Hr(NReac, NDisc) 			as heat_reaction;
	Re(NDisc+1)  				as Real					(Brief = "Reynolds Number Profile",Lower = 1E-6);
	mu(NDisc+1)  				as viscosity			(Brief = "Viscosity Profile" , Symbol = "\mu");
	fns(NDisc+1) 					as fricfactor 		(Brief = "No Slip Friction Factor");

EQUATIONS

"Inlet boundary"
	str(1).F = Inlet.F;
	str(1).T = Inlet.T;
	str(1).P = Inlet.P;
	str(1).z = Inlet.z;

"Outlet boundary"
	Outlet.F = str(NDisc+1).F;
	Outlet.T = str(NDisc+1).T;
	Outlet.P = str(NDisc+1).P;
	Outlet.z = str(NDisc+1).z;
	Outlet.h = str(NDisc+1).h;
	Outlet.v = str(NDisc+1).v;

"Reactor Initial Length"
	Lincr(1) = 0*'m';

	C(:,1)*vol(1) = Inlet.z;
	
for i in [1:NDisc] do

for c in [1:NComp] do

"Component Molar Balance"
	#diff(M(c,i)) = str(i).F*str(i).z(c) - str(i+1).F*str(i+1).z(c) + sum(stoic(c,:)*r(:, i)) * dv;
	diff(C(c,i+1)*dv) = (vel(i)*C(c,i) - vel(i+1)*C(c,i+1) )*Across + sum(stoic(c,:)*r(:, i))*dv;

end

"Constraint"
	sum(str(i+1).z ) = 1;
#	Mt(i) = sum(M(:,i));

"Molar concentration"
	#C(:,i) = rhom(i)*str(i).z;
	str(i+1).z = C(:,i+1) * vol(i+1);

#"Molar fraction"
	#str(i+1).z * Mt(i) = M(:,i);

#"Geometrical constraint"
	#Mt(i)  = dv * rhom(i);

"Energy Balance"
	diff(E(i)) = str(i).F*str(i).h - str(i+1).F*str(i+1).h +	sum(Hr(:,i)*r(:,i)) * dv - q(i);

"Energy Holdup"
	#E(i) = Mt(i)*str(i+1).h - str(i+1).P*dv;
	E(i) = dv/vol(i+1) * str(i+1).h - str(i+1).P*dv;

"Outlet Pressure"
	str(i+1).P =str(1).P - dP(i+1);

"Incremental Length"
	Lincr(i+1) = i*dx;

end

for i in [1:NDisc+1] do

"Incremental Pressure Drop"
	dP(i) = 0.5*fns(i)*Lincr(i)*rho(i)*vel(i)*vel(i)/Dpipe; # simple equation for pressure drop (Darcy Equation)

"Specific Volume"
	vol(i) = PP.VapourVolume(str(i).T, str(i).P, str(i).z);

"Specific Mass"
	rho(i) = PP.VapourDensity(str(i).T, str(i).P, str(i).z);

"Molar Density"
	rhom(i)* vol(i) = 1;

"Viscosity" 
	mu(i) = 0.027*'cP';#PP.VapourViscosity(str(i).T,str(i).P,str(i).z); # VRTherm mu=16 cP !!!!!!

"Reynolds Number"
	Re(i)*mu(i)	= rho(i)*vel(i)*Dpipe;
	
"Velocity"
	str(i).F = vel(i)*Across*rhom(i);

if Re(i) > 2300
	
	then
"Friction Factor  for Pressure Drop - Turbulent Flow"
	#1/sqrt(fns(i))= -2*log(abs(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i))));
	1/sqrt(fns(i))= -2*log(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i)));
	
	else
"Friction Factor for Pressure Drop - laminar Flow"
	fns(i)*Re(i) = 16;

end

end

end