#*-------------------------------------------------------------------
* 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 735 2009-02-27 02:00:39Z bicca $
*--------------------------------------------------------------------#

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

	DEVICES
	FEED as stream_cstr;
	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
#	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';
	
	PIDL.PID_Select = "Ideal_AWBT";
	PIDT.PID_Select = "Ideal_AWBT";
	
	EQUATIONS

	"Dimensionless level to connect PID"
	L_ad*(Lmax-Lmin)=CSTR1.h-Lmin;
	PIDL.Ports.input=L_ad;

	"Dimensionless temperature to connect PID"
	T_ad*(Tmax-Tmin)=CSTR1.T-Tmin;
	PIDT.Ports.input=T_ad;
	
	"Manipulated Variables"
	CSTR1.x  = PIDL.Ports.output;
	CSTR1.Tw  = PIDT.Ports.output*(Tmax-Tmin)+Tmin;
	
#   Level control: PID parameters
	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';

#   Temperature control: PID parameters
	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';	
	
	"Operating range for control variables"
	Lmax=5*'m';
	Lmin=0*'m';
	Tmax=700*'K';
	Tmin=230*'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';
	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