source: branches/gui/sample/miscellaneous/sample_high_index_optimal.mso @ 928

#*-------------------------------------------------------------------
* 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 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