#*-------------------------------------------------------------------
* 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 CSTR
*----------------------------------------------------------------------
* Author: Argimiro R. Secchi
* $Id: CSTR_noniso.mso 295 2008-06-15 19:48:17Z arge $
*--------------------------------------------------------------------*#

using "types";
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 uncontrolled CSTR and multiple steady-states
FlowSheet CSTR_no_control

	DEVICES
	FEED as stream_cstr;
	CSTR1 as CSTR;

	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';
	
	EQUATIONS
	"Manipulated Variables"
	CSTR1.x  = 1;
	CSTR1.Tw  = 300 * 'K';

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

#   Disturbance
	if time < 50 * 'h' then	
		"Feed Temperature" FEED.T = 300 * 'K';
	else
		"Feed Temperature" FEED.T = 350 * 'K';
	end

	INITIAL
	"Concentration" CSTR1.Ca = 50 * 'kmol/m^3';
	"Level" CSTR1.h = 1.7 * 'm';
	"Temperature" CSTR1.T = 570 * 'K'; # increase to 580 K to change steady-state
	
	OPTIONS
	TimeStep = 1;
	TimeEnd = 100;
	TimeUnit = 'h';
	DAESolver(File = "dassl");
end

Model CSTR_pid as CSTR
	
	VARIABLES
out	LI as control_signal (Brief="Level Indicator");
out	TI as control_signal (Brief="Temperature Indicator");
	
	EQUATIONS
	"Level sensor"
	LI * 'm' = h;

	"Temperature sensor"
	TI * 'K' = T; 
	
end

# Process with controlled CSTR and multiple steady-states
FlowSheet CSTR_controller

	DEVICES
	FEED as stream_cstr;
	CSTR1 as CSTR_pid;
	PIDL as PID;
	PIDT as PID;
 
	VARIABLES
	Lsp as length;
	Tsp as temperature;


	CONNECTIONS
	FEED to CSTR1.Inlet;
	CSTR1.LI to PIDL.Input;
	CSTR1.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';
	
	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';

#   Temperature control: PID parameters
	PIDT.bias = 0;
	PIDT.alpha=0.1;
	PIDT.Action="Direct";
	PIDT.gamma=1;
	PIDT.beta=1;
	PIDT.Clip="Clipped";
	PIDT.Mode="Automatic";
	PIDT.gain=1;
	PIDT.intTime=2.5*'h';
	PIDT.derivTime=1*'h';
	PIDT.tau=1*'s';
	PIDT.tauSet=1*'s';	
	
#	Operating range for control variables
	PIDL.MaxInput=5;
	PIDL.MinInput=0;
	PIDT.MaxInput=700;
	PIDT.MinInput=230;	
	PIDT.MaxOutput=700;
	PIDT.MinOutput=230;	

	EQUATIONS

	"Setpoints"
	PIDL.SetPoint * 'm' = Lsp;
	PIDT.SetPoint * 'K' = Tsp;

	"Manipulated Variables"
	CSTR1.x = PIDL.Output;
	CSTR1.Tw = PIDT.Output * '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