#*-------------------------------------------------------------------
* 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 column model
*--------------------------------------------------------------------
* 
* This sample file needs VRTherm DEMO(www.vrtech.com.br) to run
* SectionColumn_Test and needs VRTherm full to run the distillation
* column flowsheet.
*
*----------------------------------------------------------------------
* Author: Paula B. Staudt
* $Id: sample_column.mso 112 2007-01-15 01:01:46Z arge $
*--------------------------------------------------------------------*#

using "stage_separators/column";
using "controllers/PIDIncr";

# column section with 2 trays
FlowSheet SectionColumn_Test_with2tray
	PARAMETERS
	PP	as CalcObject(Brief="Physical Properties",File="vrpp");
	NComp	as Integer;

	SET
	PP.Components = [ "isobutane", "benzene"];
	PP.LiquidModel = "PR";
	PP.VapourModel = "PR";
	NComp = PP.NumberOfComponents;
	
	DEVICES
	sec as Section_Column;
	feed as stream_therm;
	reb as stream_therm;
	cond as stream_therm;
	zero as stream;
	
	CONNECTIONS
	feed to sec.trays(2).Inlet;
	zero to sec.trays(1).Inlet;
	reb to sec.trays(2).InletV;
	cond to sec.trays(1).InletL;
	
	SPECIFY
	feed.F = 113.4 * "kmol/h";
	feed.T = 291 * "K";
	feed.P = 168.3 * "kPa";
	feed.z = [0.5, 0.5];
	feed.v = 0;

	zero.F = 0 * "kmol/h";
	zero.T = 300 * "K";
	zero.P = 1 * "atm";
	zero.z = [0.5, 0.5];
	zero.v = 0;
	zero.h = 0 * "J/mol";
	
	cond.F = 68 * "kmol/h";
	cond.P = 150 * "kPa";
	cond.T = 281.75 * "K";
	cond.v = 0.0;
	cond.z = [0.6664, 0.3336];

	reb.F = 153 * "kmol/h";
	reb.P = 185 * "kPa";
	reb.T = 328.12 * "K";
	reb.z = [0.001848, 0.9982];
	reb.v = 1.0;
	
	sec.trays.Emv = 1;

	SET
	sec.NTrays = 2;
	#COLUMN
	sec.trays.V = 4 * "ft^3";
	sec.trays.Ah = 0.394 * "ft^2";
	sec.trays.lw = 20.94 * "in";
	sec.trays.hw = 0.125 * "ft";
	sec.trays.Q = 0 * "kW";
	sec.trays.beta = 0.6;
	sec.trays.alfa = 4;
	sec.trays.Ap = 3.94 * "ft^2";

	INITIAL
	sec.trays.OutletL.T = 290 *"K";
	sec.trays.Level = 0.9 * sec.trays.hw;
	sec.trays.OutletL.z(1) = 0.5;
	
	OPTIONS
	#relativeAccuracy = 1e-3;
	NLASolver = "sundials";
	time = [0:10:1000];
	
# After running few seconds of transient the steady-state
# can be obtained by using the results from that transient:

	#guessFile="SectionColumn_Test_with2tray.rlt";
	#mode = "steady";
end

# column section with 8 trays
FlowSheet SectionColumn_Test_with8tray
	PARAMETERS
	PP	as CalcObject(Brief="Physical Properties",File="vrpp");
	NComp	as Integer;

	SET
	PP.Components = [ "isobutane", "benzene"];
	PP.LiquidModel = "PR";
	PP.VapourModel = "PR";
	PP.Derivatives = 0;
	NComp = PP.NumberOfComponents;
	
	DEVICES
	sec as Section_Column;
	feed as stream_therm;
	reb as stream_therm;
	cond as stream_therm;
	zero as stream;
	
	CONNECTIONS
	feed to sec.trays(5).Inlet;
	
	zero to sec.trays([1:4]).Inlet;
	zero to sec.trays([6:8]).Inlet;
	
	reb to sec.trays(8).InletV;
	cond to sec.trays(1).InletL;
	
	SPECIFY
	feed.F = 113.4 * "kmol/h";
	feed.T = 291 * "K";
	feed.P = 168.3 * "kPa";
	feed.z = [0.5, 0.5];
	feed.v = 0;

	zero.F = 0 * "kmol/h";
	zero.T = 300 * "K";
	zero.P = 1 * "atm";
	zero.z = [0.5, 0.5];
	zero.v = 0;
	zero.h = 0 * "J/mol";
	
	cond.F = 68 * "kmol/h";
	cond.P = 150 * "kPa";
	cond.T = 281.75 * "K";
	cond.v = 0.0;
	cond.z = [0.6664, 0.3336];

	reb.F = 153 * "kmol/h";
	reb.P = 185 * "kPa";
	reb.T = 328.12 * "K";
	reb.z = [0.001848, 0.9982];
	reb.v = 1.0;
	
	sec.trays.Emv = 1;

	SET
	sec.NTrays = 8;
	#COLUMN
	sec.trays.V = 4 * "ft^3";
	sec.trays.Ah = 0.394 * "ft^2";
	sec.trays.lw = 20.94 * "in";
	sec.trays.hw = 0.125 * "ft";
	sec.trays.Q = 0 * "kW";
	sec.trays.beta = 0.6;
	sec.trays.alfa = 4;
	sec.trays.Ap = 3.94 * "ft^2";

	INITIAL
	sec.trays.OutletL.T = [290:(330-290)/(sec.NTrays-1):330] *"K";
	sec.trays.Level = 0.3 * sec.trays.hw;
	sec.trays.OutletL.z(1) = 0.5;
	
	OPTIONS
	relativeAccuracy = 1e-5;
	time = [0:1:100];
	#guessFile="SectionColumn_Test_with8tray.rlt";
	#mode = "steady";
end


FlowSheet Distillation_kettle_cond_Test
	PARAMETERS
	PP	as CalcObject(Brief="Physical Properties",File="vrpp");
	NComp	as Integer;
	
	VARIABLES
	Qc as heat_rate (Brief="Heat rate removed from condenser");
	Qr as heat_rate (Brief="Heat rate supplied to reboiler");
	
	SET
	PP.Components = [ "isobutane", "n-pentane", "propylene",
		"benzene", "isobutene" ];
	PP.LiquidModel = "PR";
	PP.VapourModel = "PR";
	PP.Derivatives = 1;
	NComp = PP.NumberOfComponents;

	DEVICES
	col as Distillation_kettle_cond;
	feed as streamTP;
	zero as stream;
	
	CONNECTIONS
	feed to col.trays(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.trays([1:4]).Inlet;
	zero to col.trays([6:col.NTrays]).Inlet;
	Qc to col.cond.Q;
	Qr to col.reb.Q;
	
	SPECIFY
	feed.F = 113.4 * "kmol/h";
	feed.T = 291 * "K";
	feed.P = 168.3 * "kPa";
	feed.z = 1/NComp;
	
	zero.F = 0 * "kmol/h";
	zero.T = 300 * "K";
	zero.P = 1 * "atm";
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * "J/mol";
	
	col.sptop.Outlet2.F = 85 * "kmol/h";
	col.reb.OutletL.F = 28.4 * "kmol/h";
	col.sptop.frac = 0.444445;
	col.cond.OutletV.F = 0 * "kmol/h";
	Qr = 3.7743e6 * "kJ/h";
	Qc = -3.71e6 * "kJ/h";
	col.pump1.dP = 16 * "kPa";
	col.trays.Emv = 1;
	
	SET
	col.NTrays = 8;
	col.cond.V = 2 * "m^3";
	col.cond.Across = 1 * "m^2";
	col.trays.V = 4 * "ft^3";
	col.trays.Ah = 0.394 * "ft^2";
	col.trays.lw = 20.94 * "in";
	col.trays.hw = 0.125 * "ft";
	col.trays.Q = 0 * "kW";
	col.trays.beta = 0.6;
	col.trays.alfa = 4;
	col.trays.Ap = 3.94 * "ft^2";
	col.reb.V = 2 * "m^3";
	col.reb.Across = 1 * "m^2";
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *"K";
	col.cond.Level = 1 * "m";
	col.cond.OutletL.z([1:4]) = [0.65, 0.05, 0.01, 0.01];

	# reboiler
	col.reb.OutletL.T = 330 *"K";
	col.reb.Level = 1 * "m";
	col.reb.OutletL.z([1:4]) = [0.1, 0.7, 0.01, 0.01];

	# column trays
	col.trays.OutletL.T = [290:(330-290)/(col.NTrays-1):330] * "K";
	col.trays.Level = 1.2 * col.trays.hw;
	col.trays.OutletL.z([1:4]) = [0.5, 0.05, 0.01, 0.01];

	OPTIONS
	relativeAccuracy = 1e-3;
	time = [0:0.01:1, 2:50];
	#guessFile="Distillation_kettle_cond_Test.rlt";
	#mode = "steady";	
end

FlowSheet Column_ctrl
	PARAMETERS
	PP	as CalcObject(Brief="Physical Properties",File="vrpp");
	NComp	as Integer;
	
	Qcmin as heat_rate (Brief="Minimum Condenser Heat supplied");
	Qcmax as heat_rate (Brief="Maximum Condenser Heat supplied");
	Qrmin as heat_rate (Brief="Minimum Reboiler Heat supplied");
	Qrmax as heat_rate (Brief="Maximum Reboiler Heat supplied");
	Frmin as flow_mol (Brief="Minimum bottom flow rate");
	Frmax as flow_mol (Brief="Maximum bottom flow rate");
	Fcmin as flow_mol (Brief="Minimum reflux flow rate");
	Fcmax as flow_mol (Brief="Maximum reflux flow rate");
	Hmint as length (Brief="Minimum liquid level in top tank");
	Hmaxt as length (Brief="Maximum liquid level in top tank");
    Hminb as length (Brief="Minimum liquid level in reboiler");
	Hmaxb as length (Brief="Maximum liquid level in reboiler");
	Pmax as pressure (Brief="Maximum column pressure");
	Pmin as pressure (Brief="Minimum column pressure");
	Tmax as temperature (Brief="Maximum column temperature");
	Tmin as temperature (Brief="Minimum column temperature");

	VARIABLES
	Qc as heat_rate (Brief="Heat rate removed from condenser");
	Qr as heat_rate (Brief="Heat rate supplied to reboiler");
	Had_top as Real (Brief="Dimensionless condenser level");
	Had_bot as Real (Brief="Dimensionless reboiler level");
	Pad as Real (Brief="Dimensionless pressure");
	Tad as Real (Brief="Dimensionless temperature");
	RR	as positive (Brief="Reflux ratio");

	SET
	PP.Components = [ "isobutane", "n-pentane", "propylene",
		"benzene", "isobutene" ];
	PP.LiquidModel = "PR";
	PP.VapourModel = "PR";
	PP.Derivatives = 1;
	NComp = PP.NumberOfComponents;

	DEVICES
	col as Distillation_kettle_cond;
	feed as streamTP;
	zero as stream;
	TCcond as PIDIncr_Ideal_AW;
	LCtop as PIDIncr_Ideal_AW;
	LCbot as PIDIncr_Ideal_AW;
	PC as PIDIncr_Ideal_AW;

	CONNECTIONS
	feed to col.trays(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.trays([1:4]).Inlet;
	zero to col.trays([6:col.NTrays]).Inlet;
	Qc to col.cond.Q;
	Qr to col.reb.Q;

	EQUATIONS
   "Temperature Controller"
	TCcond.Parameters.tau = 0*"s";	
	TCcond.Parameters.tauSet = 0*"s";	
	TCcond.Parameters.alpha = 0.3;
	TCcond.Parameters.bias = 0.5;	
	TCcond.Parameters.gamma = 1;
	TCcond.Parameters.beta = 1;
	TCcond.Options.action = 1;
	TCcond.Options.clip = 1;
	TCcond.Options.autoMan = 0;
	TCcond.Parameters.intTime = 60*"s";
	TCcond.Parameters.gain = 0.6;
	TCcond.Parameters.derivTime = 1*"s";
	TCcond.Ports.setPoint = ((15+273.15) * "K" - Tmin)/(Tmax-Tmin);
	TCcond.Ports.input = Tad;
	Tad = (col.cond.OutletL.T-Tmin)/(Tmax-Tmin);
	Qc = Qcmin+(Qcmax-Qcmin)*TCcond.Ports.output;	

	"Pressure Controller"
	PC.Parameters.tau = 0*"s";	
	PC.Parameters.tauSet = 0*"s";	
	PC.Parameters.alpha = 0.3;
	PC.Parameters.bias = 0;	
	PC.Parameters.gamma = 1;
	PC.Parameters.beta = 1;
	PC.Options.action = -1;
	PC.Options.clip = 1;
	PC.Options.autoMan = 0;
	PC.Parameters.intTime = 50*"s";
	PC.Parameters.gain = 0.5;
	PC.Parameters.derivTime = 1*"s";
	PC.Ports.setPoint = (2.0*"bar"-Pmin)/(Pmax-Pmin);
	PC.Ports.input = Pad;
	Pad = (col.cond.OutletV.P-Pmin)/(Pmax-Pmin);
	col.cond.OutletV.F = (Fcmin+(Fcmax-Fcmin)*PC.Ports.output);	
	
	"Ttop Level Controller"
	LCtop.Parameters.tau = 0*"s";	
	LCtop.Parameters.tauSet = 0*"s";	
	LCtop.Parameters.alpha = 0.3;
	LCtop.Parameters.bias = 0.5;	
	LCtop.Parameters.gamma = 1;
	LCtop.Parameters.beta = 1;
	LCtop.Options.action = -1;
	LCtop.Options.clip = 1;
	LCtop.Options.autoMan = 0;
	LCtop.Parameters.intTime = 10*"s";
	LCtop.Parameters.gain = 1;
	LCtop.Parameters.derivTime = 1*"s";
	LCtop.Ports.setPoint = (1.0 * "m" - Hmint)/(Hmaxt-Hmint);
	LCtop.Ports.input = Had_top;
	Had_top = (col.cond.Level-Hmint)/(Hmaxt-Hmint);
	col.sptop.Outlet1.F = Fcmin + (Fcmax-Fcmin) * LCtop.Ports.output;

	"Tbottom Level Controller"
	LCbot.Parameters.tau = 0*"s";	
	LCbot.Parameters.tauSet = 0*"s";	
	LCbot.Parameters.alpha = 0.3;
	LCbot.Parameters.bias = 0.5;	
	LCbot.Parameters.gamma = 1;
	LCbot.Parameters.beta = 1;
	LCbot.Options.action = -1;
	LCbot.Options.clip = 1;
	LCbot.Options.autoMan = 0;
	LCbot.Parameters.intTime = 100*"s";
	LCbot.Parameters.gain = 1;
	LCbot.Parameters.derivTime = 1*"s";
	LCbot.Ports.setPoint = (1.0 * "m" - Hminb)/(Hmaxb-Hminb);
	LCbot.Ports.input = Had_bot;
	Had_bot = (col.reb.Level-Hminb)/(Hmaxb-Hminb);
	col.reb.OutletL.F = Frmin + (Frmax-Frmin) * LCbot.Ports.output;

	RR * (col.cond.OutletV.F + col.sptop.Outlet1.F) = col.sptop.Outlet2.F;
	
	if time < 1 * "h" then
		col.sptop.Outlet2.F = 70 * "kmol/h"; # reflux
	else
		col.sptop.Outlet2.F = 85 * "kmol/h"; # reflux
	end

	SPECIFY
	feed.F = 113.4 * "kmol/h";
	feed.T = 291 * "K";
	feed.P = 168.3 * "kPa";
	feed.z = 1/NComp;
	
	zero.F = 0 * "kmol/h";
	zero.T = 300 * "K";
	zero.P = 1 * "atm";
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * "J/mol";
	
	Qr = 3e6 * "kJ/h";
	col.pump1.dP = 16 * "kPa";
	col.trays.Emv = 1;

	SET
	col.NTrays = 8;
	col.cond.V = 2 * "m^3";
	col.cond.Across = 1 * "m^2";
	col.trays.V = 4 * "ft^3";
	col.trays.Ah = 0.394 * "ft^2";
	col.trays.lw = 20.94 * "in";
	col.trays.hw = 0.125 * "ft";
	col.trays.Q = 0 * "kW";
	col.trays.beta = 0.6;
	col.trays.alfa = 4;
	col.trays.Ap = 3.94 * "ft^2";
	col.reb.V = 2 * "m^3";
	col.reb.Across = 1 * "m^2";

	Qrmax = 5e6 * "kJ/h";
	Qrmin = 1e6 * "kJ/h";
	Frmin = 0 * "kmol/h";
	Frmax = 60 * "kmol/h";
	Fcmin = 0 * "kmol/h";
	Fcmax = 120 * "kmol/h";
	Hmint = 0 * "m";
	Hmaxt = 2 * "m";
	Hminb = 0 * "m";
	Hmaxb = 2 * "m";
	Pmin = 0.5 * "bar";
	Pmax = 4 * "bar";
	Qcmax = -5e5 * "kJ/h";
	Qcmin = -5e6 * "kJ/h";
	Tmax = (30+273.15) * "K";
	Tmin = (-20+273.15) * "K";
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *"K";
	col.cond.Level = 1 * "m";
	col.cond.OutletL.z([1:4]) = [0.2, 0.2, 0.4, 0.05];

	# reboiler
	col.reb.OutletL.T = 350 *"K";
	col.reb.Level = 1 * "m";
	col.reb.OutletL.z([1:4]) = [0.1, 0.4, 0.1, 0.3];
	
	# column trays
	col.trays.OutletL.T = [290:(330-290)/(col.NTrays-1):330] * "K";
	col.trays.Level = 1.2 * col.trays.hw;
	col.trays.OutletL.z([1:4]) = [0.15, 0.3, 0.25, 0.2];

	OPTIONS
	relativeAccuracy = 1e-3;
	time = [0:0.01:0.1, 0.11:0.01:2]*"h";
	#initialFile = "Column_ctrl.rlt";
	#guessFile = "Column_ctrl.rlt";
	#mode = "steady";
end