using "controllers/PIDS";

const_valv as positive(Brief = "Valve Constant", Default=1,Lower=0,Upper=100, Unit="m^2.5/h");

Model corrente
	VARIABLES
	Ca     as conc_mol;
	F      as flow_vol;
	T      as temperature;
end

Model CSTR

	PARAMETERS
	ko  as frequency		(Unit="1/h");
	A   as area;
	At  as area;	
	Ea  as energy_mol 		(Unit="kJ/kmol");
	R   as Real 			(Unit="kJ/mol/K");
	ro  as dens_mol 		(Unit="kmol/m^3");
	Cp  as cp_mol 			(Unit="kJ/kmol/K");
	U   as heat_trans_coeff	(Unit="kW/m^2/K");
	Hr  as heat_reaction    (Unit="kJ/kmol");
	
	VARIABLES
	Cv       as const_valv;
	T        as temperature;
	Tw       as temperature;
	V        as volume;
	Ca       as conc_mol;
	h        as length;
	tau      as time_h;
	rA	     as reaction_mol;
	k        as frequency	(Unit="1/h"); 
in  Inlet    as corrente;
out Outlet   as corrente;

	EQUATIONS

	"Balanço de Massa Global"
	diff(V) = Inlet.F - Outlet.F;
	"Balanço de Massa por Componente"
	tau * diff(Ca) = (Inlet.Ca - Ca) - (-rA) * tau;
	"Mistura perfeita"
	Outlet.Ca = Ca;
	Outlet.T  = T;
	"Taxa de reação"
	-rA = k * Ca;
	"Equação de Arrhenius"
	k=ko*exp(-Ea/(R*T));
	"Tempo de residência médio"
	tau * Inlet.F = V;
	"Geometria"
	A * h = V;
	"Equação da válvula"
	Outlet.F = Cv * sqrt(h);
	"Balanço de energia"
	tau * diff(T)=(Inlet.T - T) - U*At*(T-Tw)/(ro*V*Cp)*tau + (-Hr)*(-rA)*tau/(ro*Cp);

end 

# processo com 1 CSTR controlado
FlowSheet CSTR_controller

	DEVICES
	FEED as corrente;
	CSTR as CSTR;
	PIDL as PID_Ideal_AWBT;
	PIDT as PID_Ideal_AWBT;
 
	VARIABLES
	L_ad as Real;
	Lmin as length;
	Lmax as length;
	T_ad as Real;
	Tmin as temperature;
	Tmax as temperature;
	
	CONNECTIONS
	FEED to CSTR.Inlet;
	
	SET
#	Parâmetros do CSTR"
	CSTR.R   = 8.3144 * "kJ/kmol/K";
	CSTR.U   = 300 * "kJ/h/m^2/K";
	CSTR.ro  = 55.56 * "kmol/m^3";
	CSTR.Cp  = 70*"kJ/kmol/K";
	CSTR.Hr  = -7000 * "kJ/kmol";
	CSTR.Ea  = 6e4 * "kJ/kmol";
	CSTR.ko  = 89 * "1/s";
	CSTR.A   = 8 * "m^2";
	CSTR.At  = 25 * "m^2";
	
	EQUATIONS

	"Equações do PID para controle de nível"
	L_ad*(Lmax-Lmin)=CSTR.h-Lmin;
	PIDL.Ports.input=L_ad;

	"Equações do PID para controle de temperatura"
	T_ad*(Tmax-Tmin)=CSTR.T-Tmin;
	PIDT.Ports.input=T_ad;
	
	"Variáveis manipulada"
	CSTR.Cv  = 2.2136 * "m^2.5/h"*PIDL.Ports.output;
	CSTR.Tw  = PIDT.Ports.output*(Tmax-Tmin)+Tmin;
	
	#distúrbio regulatório
	if time<1.6e5 then	
		FEED.T = 300 * "K";
	else
		FEED.T = 285 * "K";
	end


	#Parâmetros do PID de nível
	PIDL.Parameters.bias = 0;
	PIDL.Parameters.alpha=0.1;
	PIDL.Options.action=-1;
	PIDL.Parameters.gamma=1;
	PIDL.Parameters.beta=1;
	PIDL.Options.clip=1;
	PIDL.Options.autoMan=0;
	PIDL.Parameters.gain=0.15;
	PIDL.Parameters.intTime=100*"s";
	PIDL.Parameters.derivTime=1*"s";
	PIDL.Ports.setPoint=0.55;
	PIDL.Parameters.tau=1*"s";
	PIDL.Parameters.tauSet=1*"s";
	
	PIDT.Parameters.bias = 0;
	PIDT.Parameters.alpha=0.1;
	PIDT.Options.action=-1;
	PIDT.Parameters.gamma=1;
	PIDT.Parameters.beta=1;
	PIDT.Options.clip=1;
	PIDT.Options.autoMan=0;
	PIDT.Parameters.gain=0.85;
	PIDT.Parameters.intTime=100*"s";
	PIDT.Parameters.derivTime=1*"s";
	PIDT.Ports.setPoint=0.85;
	PIDT.Parameters.tau=1*"s";
	PIDT.Parameters.tauSet=1*"s";	
	
	"Valores limites para normalizações"
	Lmax=5*"m";
	Lmin=0*"m";
	Tmax=700*"K";
	Tmin=230*"K";	
	
	"Variáveís da corrente de alimentação"
	FEED.Ca = 300 * "kmol/m^3";
	FEED.F = 3.5 * "m^3/h";

	INITIAL
	CSTR.Ca = 50 * "kmol/m^3";
	CSTR.h = 2.5 * "m";
	CSTR.T = 650 * "K";
	
	OPTIONS
	time = [0:0.1:1 1:1:100] * "h";
	outputLevel = "high";
end