#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2008 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 - 2008 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: Argimiro R. Secchi
* $Id: CSTR_noniso_pid.mso 536 2008-06-15 23:18:35Z arge $
*--------------------------------------------------------------------*#

using "controllers/PIDs";

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

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

Model CSTR

	PARAMETERS
	ko  as frequency		(DisplayUnit='1/h');
	D	as length;
	A   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');
	pi  as Real				(Default = 3.141593);
	Cv  as const_valv;

	VARIABLES
	At		 as area;	
	T        as temperature;
	Tw       as temperature;
	x		 as fraction;
	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 stream_cstr;
out Outlet   as stream_cstr;

	SET
	A = pi * D^2 / 4;
	
	EQUATIONS

	"Overall Mass Balance"
	diff(V) = Inlet.F - Outlet.F;
	
	"Component Mass Balance"
	V * diff(Ca) = Inlet.F * (Inlet.Ca - Ca) - (-rA) * V;

	"Average Residence Time"
	tau * Inlet.F = V;

	"Energy Balance"
	ro * V * Cp * diff(T) = Inlet.F * ro * Cp * (Inlet.T - T) + qr - q;

	"Heat Transfer Rate"
	q = U * At * (T - Tw);

	"Reaction Heat Rate"
	qr = (-Hr) * (-rA) * V;
	
	"Reaction Rate"
	-rA = k * Ca;
	
	"Arrhenius Equation"
	k = ko * exp(-Ea/(R*T));
	
	"Geometry"
	A * h = V;
	At = A + pi*D*h;
	
	"Valve Equation"
	Outlet.F = x * Cv * sqrt(h);

	"Perfect Mixture"
	Outlet.Ca = Ca;
	Outlet.T  = T;
end

# Process with controlled CSTR and multiple steady-states
FlowSheet CSTR_controller
	PARAMETERS
	Lmin as length;
	Lmax as length;
	Tmin as temperature;
	Tmax as temperature;

	DEVICES
	FEED as stream_cstr;
	CSTR1 as CSTR;
	PIDL as PID;
	PIDT as PID;
 
	VARIABLES
	Lsp as length;
	Tsp as temperature;
out	TI	as control_signal;
out	LI	as control_signal;
	
	CONNECTIONS
	FEED to CSTR1.Inlet;
	LI	 to PIDL.Input;
	TI	 to PIDT.Input;

	SET
#	CSTR Parameters
	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.D   = 3.2 * 'm';
	CSTR1.Cv  = 2.7 * 'm^2.5/h';

#	Operating range for control variables
	Lmax=5*'m';
	Lmin=0*'m';
	Tmax=700*'K';
	Tmin=230*'K';

	PIDL.PID_Select = "Ideal_AWBT";
	PIDT.PID_Select = "Ideal_AWBT";

#   Level control: PID parameters
	PIDL.bias=0;
	PIDL.alpha=0.1;
	PIDL.Action="Direct";
	PIDL.gamma=1;
	PIDL.beta=1;
	PIDL.Clip="Clipped";
	PIDL.Mode="Automatic";
	PIDL.gain=1;
	PIDL.intTime=2.5*'h';
	PIDL.derivTime=0*'s';
	PIDL.tau=1*'s';
	PIDL.tauSet=1*'s';
	PIDL.MinInput=Lmin/'m';
	PIDL.MaxInput=Lmax/'m';

#   Temperature control: PID parameters
	PIDT.bias = 0;
	PIDT.alpha=0.1;
	PIDT.Action="Reverse";
	PIDT.gamma=1;
	PIDT.beta=1;
	PIDT.Clip="Clipped";
	PIDT.Mode="Automatic";
	PIDT.gain=1;
#	PIDT.gain=15; # for a setpoint of 523K or a faster response
	PIDT.intTime=2.5*'h';
	PIDT.derivTime=1*'h';
	PIDT.tau=1*'s';
	PIDT.tauSet=1*'s';
	PIDT.MinInput=Tmin/'K';
	PIDT.MaxInput=Tmax/'K';
	PIDT.MinOutput=Tmin/'K';
	PIDT.MaxOutput=Tmax/'K';
	
	EQUATIONS

	"Dimensionless level to connect PID"
	LI=CSTR1.h/'m';

	"Dimensionless temperature to connect PID"
	TI=CSTR1.T/'K';

	"Manipulated Variables"
	CSTR1.x  = PIDL.Output;
	CSTR1.Tw  = PIDT.Output*'K';

	"Level setpoint"
	PIDL.SetPoint=Lsp/'m';

	"Temperature setpoint"
	PIDT.SetPoint=Tsp/'K';

	"Feed Stream"
	FEED.Ca = 300 * 'kmol/m^3';
	FEED.F = 3.5 * 'm^3/h';

#   Disturbance
	if time < 50 * 'h' then	
		FEED.T = 300 * 'K';
	else
		FEED.T = 285 * 'K'; # change to 350 K to saturate controller
	end

#   Set-point changes
	if time < 100 * 'h' then
	  Tsp = 630 * 'K';
	else
	  Tsp = 400 * 'K';
#	  Tsp = 523 * 'K'; # near unstable steady-state (change also the controller's gain)
	end

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

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