#*---------------------------------------------------------------------
* This file is property of the author and cannot be used, copyed
* or modified without permission.
*
* Copyright (C) 2004  the author
*----------------------------------------------------------------------
* Author: Rafael de P. Soares
*         Paula B. Staudt
* $Id: cstr.mso 1 2006-06-20 17:33:53Z rafael $
*--------------------------------------------------------------------*#

using "streams";

#*
* Generic PFR model with constant mass holdup
*
* Requires the information of:
* - Reaction values
* - Heat of reaction
* - Pressure profile
*#
Model pfr
	PARAMETERS
ext PP   as CalcObject (Brief = "External Physical Properties");
ext	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");

	SET
	Mw   = PP.MolecularWeight();

	VARIABLES
in	Inlet  as stream;
out	Outlet as stream;
	str(NDisc+1) as stream_therm;
	vel(NDisc+1) as velocity;
	vol(NDisc+1) as vol_mol;
	rho(NDisc+1) as dens_mass;
	
	q(NDisc)    as heat_rate;
	M(NComp, NDisc)    as mol (Brief = "Molar holdup");
	C(NComp, NDisc)  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;

	EQUATIONS
	"Vapourisation Fraction"
	str.v = Inlet.v;

	"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;
	
	for z in [1:NDisc]
		for c in [1:NComp]
			"Component Molar Balance"
			diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
				+ sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
		end

		"Energy Balance"
		diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
			sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);

		"Energy Holdup"
		E(z) = sum(M(:,z))*str(z+1).h - str(z+1).P*Across*L/NDisc;

		"mass flow is considered constant"
		str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar)
#		rho(z+1)*vel(z+1) = rho(z)*vel(z);  # FIXME: this is correct! But does not converge.
	
		"Molar concentration"
		C(:,z) * Across*L/NDisc = M(:,z);

		"Molar fraction"
		str(z+1).z = C(:,z) * vol(z+1);
	end	

	for z in [1:NDisc+1]
		"Specific Volume"
		vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z);
		
		"Specific Mass"
		rho(z) = sum(str(z).z*Mw)/vol(z);

		"Velocity"
		vel(z)*Across = str(z).F*vol(z);
	end
end