source: trunk/eml/reactors/gibbs.mso @ 417

#*---------------------------------------------------------------------
* 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 a Gibbs reactor
*----------------------------------------------------------------------
*
*   Description:

*       Thermodynamic equilibrium modeling of a reactor using Gibbs
*       free energy minimization approach.

*

*   Assumptions:

*       * thermodynamic equilibrium

*               * steady-state

*

*       Specify:

*               * inlet stream
*               * number of elements related to components

*               * matrix of elements by components

*               * equilibrium temperature
*
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* $Id$
*--------------------------------------------------------------------*#

using "tank_basic";


#*---------------------------------------------------------------------
*       only vapour phase
*--------------------------------------------------------------------*#
Model gibbs_vap as tank_vap
        ATTRIBUTES
        Pallete = true;
        Icon    = "icon/cstr"; 
        Brief   = "Model of a generic vapour-phase Gibbs CSTR";
        Info    = "
Requires the information of:
* number of elements
* matrix of elements (elements by compoments)
";

        PARAMETERS
outer NElem                     as Integer              (Brief="Number of elements", Default=1);
        Rg                              as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
        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; # Outlet stream

        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier");
        activ(NComp)    as Real                 (Brief="Activity", Lower=1e-20);
        phi(NComp)              as fugacity             (Brief="Fugacity coefficient", Default=1);
        
        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Default=0);
        Fi(NComp)               as flow_mol             (Brief="Component molar flow rate");

        EQUATIONS
        "Outlet stream"
        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
        
        "Mechanical equilibrium"
        Outlet.P = Outletm.P;

        "Steady-state"
        Outlet.F = sum(Fi);

        "Component molar flow rate"
        Fi = Outlet.F*Outlet.z;
        
        "Energy balance"
        Outlet.F*Outlet.h = Outletm.F*Outletm.h;

        "Element balance"
        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);

        "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 i in [1:NComp]
        "Lagrangian multiplier"
                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));

          if (Outletm.z(i) > 0) then
            "Molar conversion"
            Fi(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
        
        "Fugacity coefficient"
        phi = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
        
        "Activity"
        activ = phi*Outlet.P*Outlet.z/fs;
end


#*---------------------------------------------------------------------
*       only liquid phase
*--------------------------------------------------------------------*#
Model gibbs_liq as tank_liq
        ATTRIBUTES
        Pallete = true;
        Icon    = "icon/cstr"; 
        Brief   = "Model of a generic liquid-phase Gibbs CSTR";
        Info    = "
Requires the information of:
* number of elements
* matrix of elements (elements by compoments)
";

        PARAMETERS
outer NElem                     as Integer              (Brief="Number of elements", Default=1);
        Rg                              as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
        Ps                              as pressure             (Brief="Pressure of standard state", Default=1, DisplayUnit='atm');
        To                              as temperature  (Brief="Reference temperature", Default=298.15);
        
        VARIABLES
out Outlet                      as liquid_stream; # Outlet stream

        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier");
        activ(NComp)    as Real                 (Brief="Activity", Lower=0);
        gamma(NComp)    as fugacity             (Brief="Activity coefficient", Default=1);

        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Default=0);
        Fi(NComp)               as flow_mol             (Brief="Component molar flow rate");

        EQUATIONS
        "Outlet stream"
        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
        
        "Mechanical equilibrium"
        Outlet.P = Outletm.P;

        "Steady-state"
        Outlet.F = sum(Fi);

        "Component molar flow rate"
        Fi = Outlet.F*Outlet.z;
        
        "Energy balance"
        Outlet.F*Outlet.h = Outletm.F*Outletm.h;

        "Element balance"
        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
        
        "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 i in [1:NComp]
        "Lagrangian multiplier"
                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));

          if (Outletm.z(i) > 0) then
            "Molar conversion"
            Fi(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 coefficient"
        gamma = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
        
        "Activity"
        activ = gamma*Outlet.z*exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
end