#*-------------------------------------------------------------------
* 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.
*
*--------------------------------------------------------------------
* Sample file for controllers
*----------------------------------------------------------------------
* Author:
* $Id: CSTR_noniso_pid.mso 295 2007-06-21 19:48:17Z bicca $
*--------------------------------------------------------------------*#

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		(DisplayUnit='1/h');
	A   as area;
	At  as area;	
	Ea  as energy_mol 		(DisplayUnit='kJ/kmol');
	R   as Real 			(Unit='kJ/mol/K');
	ro  as dens_mass 		(DisplayUnit='kg/m^3');
	Cp  as cp_mass 			(DisplayUnit='kJ/kg/K');
	U   as heat_trans_coeff	(DisplayUnit='kW/m^2/K');
	Hr  as heat_reaction    (DisplayUnit='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	(DisplayUnit='1/h');
	q		 as heat_rate	(DisplayUnit='kJ/h');
	qr		 as heat_rate	(DisplayUnit='kJ/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"
	V * diff(Ca) = Inlet.F * (Inlet.Ca - Ca) - (-rA) * V;

	"Tempo de residência médio"
	tau * Inlet.F = V;

	"Balanço de energia"
	ro * V * Cp * diff(T) = Inlet.F * ro * Cp * (Inlet.T - T) + qr - q;

	"Taxa de calor transferido"
	q = U * At * (T - Tw);

	"Taxa de calor gerado"
	qr = (-Hr) * (-rA) * V;
	
	"Taxa de reação"
	-rA = k * Ca;
	
	"Equação de Arrhenius"
	k = ko * exp(-Ea/(R*T));
	
	"Geometria"
	A * h = V;
	
	"Equação da válvula"
	Outlet.F = Cv * sqrt(h);

	"Mistura perfeita"
	Outlet.Ca = Ca;
	Outlet.T  = T;
end 

# processo com 1 CSTR controlado
FlowSheet CSTR_controller

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

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

	"Equações do PID para controle de temperatura"
	T_ad*(Tmax-Tmin)=CSTR1.T-Tmin;
	PIDT.Ports.input=T_ad;
	
	"Variáveis manipulada"
	CSTR1.Cv  = 2.2136 * 'm^2.5/h' * PIDL.Ports.output;
	CSTR1.Tw  = PIDT.Ports.output*(Tmax-Tmin)+Tmin;
	
	#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=1;
	PIDL.Parameters.intTime=2.5*'h';
	PIDL.Parameters.derivTime=0*'s';
	PIDL.Ports.setPoint=(Lsp - Lmin)/(Lmax - Lmin);
	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=1;
	PIDT.Parameters.intTime=2.5*'h';
	PIDT.Parameters.derivTime=1*'h';
	PIDT.Ports.setPoint=(Tsp - Tmin)/(Tmax - Tmin);
	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';

	#distúrbio regulatório
	if time < 50 * 'h' then	
		FEED.T = 300 * 'K';
	else
		FEED.T = 285 * 'K';
	end

	#mudança de set-point
	if time < 100 * 'h' then
	  Tsp = 630 * 'K';
	else
	  Tsp = 400 * 'K';
	end

	if time < 150 * 'h' then
	  Lsp = 2.5 * 'm';
	else
	  Lsp = 4 * 'm';
	end

	INITIAL
	CSTR1.Ca = 50 * 'kmol/m^3';
	CSTR1.h = 2.5 * 'm';
	CSTR1.T = 550 * 'K';
	
	OPTIONS
	TimeStep = 1;
	TimeEnd = 250;
	TimeUnit = 'h';
	DAESolver(File = "dasslc");
end