source: trunk/eml/reactors/equil.mso @ 659

#*---------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* Model of an equilibrium reactor
*----------------------------------------------------------------------
*
*   Description:
*       Thermodynamic equilibrium modeling of a reactor based on 
*       equilibrium constants approach.
*
*   Assumptions:
*               * single-phases involved
*       * thermodynamic equilibrium
*               * steady-state
*
*       Specify:
*               * inlet stream
*               * stoichiometric matrix
*               * equilibrium temperature
*
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* $Id$
*--------------------------------------------------------------------*#

using "tank_basic";


#*---------------------------------------------------------------------
*       only vapour phase
*--------------------------------------------------------------------*#
Model equil_vap as tank_vap
        ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/cstr"; 
        Brief           = "Model of a generic vapour-phase equilibrium CSTR";
        Info            = "
== Assumptions ==
* only vapour-phase
* thermodynamic equilibrium
* steady-state

== Specify ==
* inlet stream
* stoichiometric matrix
* equilibrium temperature
";

        PARAMETERS
        NReac           as Integer              (Brief="Number of reactions", Default=1);
    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix", Symbol="\nu");
        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
        fs(NComp)       as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
        To                      as temperature  (Brief="Reference temperature", Default=298.15);

        VARIABLES
out Outlet              as vapour_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");

        G(NComp)        as energy_mol   (Brief="Gibbs free-energy of formation");
        K(NReac)        as Real                 (Brief="Equillibrium constant", Lower=0, Default=1.5);
        activ(NComp)as Real             (Brief="Activity", Symbol="\hat{a}", Lower=0, Default=0.2);
        
        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
        extent(NReac) as flow_mol       (Brief="Extent of reaction", Symbol="\xi");
        conv(NComp) as Real             (Brief="Fractional conversion of component", Symbol="X", Default=0); # Lower=-1e3, Upper=1e3);
        
        EQUATIONS
        "Outlet stream"
        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
        
        "Mechanical equilibrium"
        Outlet.P = Outletm.P;
        
        "Energy balance"
        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
        
        "Steady-state"
        Outlet.F = Inlet.F + sum(sumt(stoic*extent));

        "Gibbs free-energy of formation"
        G = PP.IdealGasGibbsOfFormation(Outlet.T);

#       "Gibbs free-energy of formation without Cp correction"
#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);

        for j in [1:NReac]
        "Gibbs free energy of reaction"
                sumt(G*stoic(:,j)) = -Rg*Outlet.T*ln(K(j));
#               K(j) = exp(-sumt(G*stoic(:,j))/(Rg*Outlet.T));

        "Equilibrium constant"
                K(j) = prod(activ^stoic(:,j));
        end

        for i in [1:NComp]
        "Outlet molar fraction"
                Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
        end     

        for i in [1:NComp]
          if (Outletm.z(i) > 1e-16) then
            "Molar conversion"
            Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
          else if (Outlet.z(i) > 0) then
                        "Molar conversion"
                                conv(i) = 1;    # ?
                        else
                        "Molar conversion"
                                conv(i) = 0;    # ?
                        end
          end
        end

        "Activity"
        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
end


#*---------------------------------------------------------------------
*       only liquid-phase
*--------------------------------------------------------------------*#
Model equil_liq as tank_liq
        ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/cstr"; 
        Brief           = "Model of a generic liquid-phase equilibrium CSTR";
        Info            = "
== Assumptions ==
* only liquid-phase
* thermodynamic equilibrium
* steady-state

== Specify ==
* inlet stream
* stoichiometric matrix
* equilibrium temperature
";

        PARAMETERS
        NReac           as Integer              (Brief="Number of reactions", Default=1);
    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix", Symbol="\nu");
        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
        Ps                      as pressure     (Brief="Standard pressure", Default=1, DisplayUnit='bar');
        To                      as temperature  (Brief="Reference temperature", Default=298.15);

        VARIABLES
out Outlet              as liquid_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");

        G(NReac)        as enth_mol     (Brief="Gibbs free-energy of formation");
        K(NReac)        as fraction             (Brief="Equillibrium constant");
        activ(NComp)as Real             (Brief="Activity", Symbol="\hat{a}");
        
        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
        extent(NReac)as flow_mol        (Brief="Extent of reaction", Symbol="\xi");
        conv(NComp) as Real             (Brief="Fractional conversion of component", Symbol="X", Default=0);
        
        EQUATIONS
        "Outlet stream"
        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
        
        "Mechanical equilibrium"
        Outlet.P = Outletm.P;
        
        "Energy balance"
        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
        
        "Steady-state"
        Outlet.F = Inlet.F + sum(sumt(stoic*extent));
        
        "Gibbs free-energy of formation"
        G = PP.IdealGasGibbsOfFormation(Outlet.T);

#       "Gibbs free-energy of formation without Cp correction"
#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);

        "Gibbs free energy of reaction"
        sumt(G*stoic) = -Rg*Outlet.T*ln(K);
#       K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
        
        for j in [1:NReac]
        "Equilibrium constant"
                K(j) = prod(activ^stoic(:,j));
        end
        
        for i in [1:NComp]
        "Outlet molar fraction"
                Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
        end     

        for i in [1:NComp]
          if (Outletm.z(i) > 1e-16) then
            "Molar conversion"
            Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
          else if (Outlet.z(i) > 0) then
                        "Molar conversion"
                                conv(i) = 1;    # ?
                        else
                        "Molar conversion"
                                conv(i) = 0;    # ?
                        end
          end
        end
        
        "Activity"
        activ = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z
                *exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
end