#*-------------------------------------------------------------------
* 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 for a high-index optimal control problem.
*--------------------------------------------------------------------
* Author: Rafael de Pelegrini Soares
* $Id$
*--------------------------------------------------------------------*# 
using "types";

Model FlashRaoult
	ATTRIBUTES
	Info = "
	This is a very simple (wrong) model with dynamics only on the energy.
	
	It should be used for ilustration purposes only.";
	
	PARAMETERS
	NComp as Integer;
	
	# Antoine constants
	A(NComp) as Real;
	B(NComp) as Real;
	C(NComp) as Real;
	
	Cv as Real(Unit = 'J/mol/K');
	DHvap(NComp) as energy_mol;

	VARIABLES
	F as flow_mol;
	L as flow_mol;
	V as flow_mol;
	z(NComp) as fraction;
	x(NComp) as fraction;
	y(NComp) as fraction;
	n(NComp) as mol;
	nt       as mol;
	
	T as temperature;
	P as pressure;
	Psat(NComp) as pressure;
	
	Q as power;
	E as energy;
	
	EQUATIONS	
	"Component Molar Balance"
	diff(n) = F*z - (L*x + V*y);
	
	nt = sum(n);
	x = n/nt;
	
	"Energy Balance"
	diff(E) = Q - V*sum(DHvap*y);
	
	"Internal energy"
	E = nt*Cv*T;
	
	"Raoult's Law"
	P*y = Psat*x;
	
	"Antoine for Vapour Pressure"
	ln(Psat/'kPa') = A - B/(T/'K'- 273.15 + C);
	
	"Molar Fraction sum"
	sum(y) = 1;
end


FlowSheet FlashRaoultTest
	DEVICES
	fl as FlashRaoult;
	
	SET
	fl.NComp = 3;
	# Antoine constants (Acetone, Acetonitrile, Nitromethane)
	fl.A = [14.31, 14.89, 14.75];
	fl.B = [2756,  3413,  3331];
	fl.C = [228,   250,   227];
	
	fl.Cv = 30 * 'J/mol/K';
	fl.DHvap = [10, 20, 30] * 'kJ/mol';
	
	EQUATIONS
	# Disturb
	#if time < 0.5 * 'h' then
	#	fl.z = [0.45, 0.35, 0.2];
	#else
	#	fl.z = [0.55, 0.25, 0.2];
	#end

	SPECIFY
	# Feed condition
	fl.F = 1 * 'kmol/h';
	# Steady-state feed
	#fl.z = [0.45, 0.35, 0.2];
	# Disturb on feed composition
	fl.z = [0.55, 0.25, 0.2];

	# Desired production of y1 (index 2 - will determine the "perfect" heat profile)
	fl.V = 0.1 * 'kmol/h';
	fl.y(1) = 0.7;
	fl.nt = 1 * 'kmol';

	# Default specification (index 1 - heat is given)
	#fl.V = 0.1 * 'kmol/h';
	#fl.Q = 0.37 * 'kW';
	#fl.nt = 1 * 'kmol';
	
	# Fixed Temperature (index 2 - will determine the "perfect" heat profile)
	#fl.T = (90+273.15) * 'K';
	#fl.V = 0.1 * 'kmol/h';
	#fl.nt = 1 * 'kmol';

	# Fixed Pressure (index 2 - will determine the "perfect" heat profile)
	#fl.P = 1.2 * 'atm';
	#fl.V = 0.1 * 'kmol/h';
	#fl.nt = 1 * 'kmol';

	INITIAL
	#fl.T = (80+273.15) * 'K';
	#fl.nt = 1 * 'kmol';

	# Steady-state composition with feed = [0.45, 0.35, 0.2];
	fl.x(1) = 0.4222;
	fl.x(2) = 0.3622;

	# steady state compositions
	# diff(fl.n(1:2)) = 0 * 'mol/s';
	
	OPTIONS
	TimeEnd = 4;
	TimeStep = 0.05;
	TimeUnit = 'h';

	#DAESolver(File="dasslc"); # slow integration
	DAESolver(File="mebdf");  # much faster
	
	Dynamic = true;
end