source: trunk/sample/miscellaneous/sample_flashPH.mso @ 304

#*-------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* Author: Rafael de P. Soares and Paula B. Staudt
* $Id: sample_flashPH.mso 289 2007-06-20 15:07:08Z rafael $
*--------------------------------------------------------------------*#

using "streams";

Model FlashPHSteady
        ATTRIBUTES
        Pallete         = true;
        Icon            = "Flash"; 
        Brief           = "Model of a static PH flash.";
        Info            = "
        This model shows how to model a pressure enthalpy flash
        directly with the EMSO modeling language.

        This model is for demonstration purposes only, the flashPH
        routine available on VRTherm is much more robust.

        Assumptions:
         * perfect mixing of both phases;

        Specify: 
         * the feed stream;
         * the heat duty;
         * the outlet pressure.
        ";      
        
        PARAMETERS
outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
outer NComp as Integer;
        B as Real(Default=1000, Brief="Regularization Factor");

        VARIABLES
in      Inlet as stream (Brief="Feed Stream");
out     OutletL as liquid_stream (Brief="Liquid outlet stream");
out     OutletV as vapour_stream (Brief="Vapour outlet stream");
in      Q as heat_rate (Brief="Rate of heat supply"); 
        vfrac as fraction(Brief="Real vaporization fraction");
        vsat as Real(Lower=-0.1, Upper=1.1, Brief="Vaporization fraction if saturated");
        Tsat as temperature(Lower=173, Upper=1473, Brief="Temperature if saturated");
        xsat(NComp) as Real(Lower=0, Upper=1, Brief="Liquid composition if saturated");
        ysat(NComp) as Real(Lower=0, Upper=1, Brief="Vapour composition if saturated");
        
        zero_one as fraction(Brief="Regularization Variable");
        one_zero as fraction(Brief="Regularization Variable");

        EQUATIONS
        "Chemical equilibrium"
        PP.LiquidFugacityCoefficient(Tsat, OutletL.P, xsat)*xsat = 
                PP.VapourFugacityCoefficient(Tsat, OutletV.P, ysat)*ysat;

        "Global Molar Balance"
        Inlet.F = OutletV.F + OutletL.F;
        OutletV.F = Inlet.F * vfrac;

        "Component Molar Balance"
        Inlet.F*Inlet.z = OutletL.F*xsat + OutletV.F*ysat;
        sum(xsat) = sum(ysat);

        "Energy Balance if saturated"
        Inlet.F*Inlet.h  + Q =
                Inlet.F*(1-vsat)*PP.LiquidEnthalpy(Tsat, OutletL.P, xsat) +
                Inlet.F*vsat*PP.VapourEnthalpy(Tsat, OutletV.P, ysat);

        "Real Energy Balance"
        Inlet.F*Inlet.h  + Q =
                Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h;

        "Thermal Equilibrium"
        OutletV.T = OutletL.T;
        
        "Mechanical Equilibrium"
        OutletV.P = OutletL.P;
        
        # regularization functions
        zero_one = (1 + tanh(B * vsat))/2;
        one_zero = (1 - tanh(B * (vsat - 1)))/2;
        
        vfrac = zero_one * one_zero * vsat + 1 - one_zero;
        OutletL.z = zero_one*one_zero*xsat + (1-zero_one*one_zero)*Inlet.z;
        OutletV.z = zero_one*one_zero*ysat + (1-zero_one*one_zero)*Inlet.z;
end

FlowSheet FlashPHTest
        PARAMETERS
        PP      as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["methane", "isobutene", "benzene"],
                LiquidModel = "IdealLiquid",
                VapourModel = "Ideal"
        );
        NComp   as Integer;

        VARIABLES
        Q       as heat_rate (Brief="Heat supplied");
        
        SET
        NComp = PP.NumberOfComponents;
        
        DEVICES
        fl as FlashPHSteady;
        s1 as source;
        
        CONNECTIONS
        s1.Outlet to fl.Inlet;
        Q to fl.Q;
        
        SPECIFY
        s1.Outlet.F = 496.3 * 'kmol/h';
        s1.Outlet.T = 338 * 'K';
        s1.Outlet.P = 507.1 * 'kPa';
        s1.Outlet.z = [0.1, 0.7,0.2];
        
        fl.OutletL.P = 2.5 * 'atm';

        Q = 0 * 'kJ/h';
        #fl.OutletL.T = 315.06 * 'K';
        
        OPTIONS
        Dynamic = false;
end