Changeset 414 for trunk/eml


Ignore:
Timestamp:
Nov 25, 2007, 2:15:43 PM (15 years ago)
Author:
Rodolfo Rodrigues
Message:

Added liquid-phase calculation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/eml/reactors/gibbs.mso

    r403 r414  
    1717*----------------------------------------------------------------------
    1818*
    19 *
     19*   Description:
     20
     21*       Thermodynamic equilibrium modeling of a reactor using Gibbs
     22*       free energy minimization approach.
     23
     24*
     25
     26*   Assumptions:
     27
     28*       * thermodynamic equilibrium
     29
     30*               * steady-state
     31
     32*
     33
     34*       Specify:
     35
     36*               * inlet stream
     37*               * number of elements related to components
     38
     39*               * matrix of elements by components
     40
     41*               * equilibrium temperature
    2042*
    2143*----------------------------------------------------------------------
     
    3052*       only vapour phase
    3153*--------------------------------------------------------------------*#
    32 
    3354Model gibbs_vap as tank_vap
    3455        ATTRIBUTES
    35         Pallete         = true;
    36         Icon            = "icon/cstr";
    37         Brief           = "Model of a generic vapour-phase Gibbs CSTR";
    38         Info            = "
     56        Pallete = true;
     57        Icon    = "icon/cstr";
     58        Brief   = "Model of a generic vapour-phase Gibbs CSTR";
     59        Info    = "
    3960Requires the information of:
    4061* number of elements
    41 * matrix of elements (elements, compounds)
     62* matrix of elements (elements by compoments)
    4263";
    4364
     
    4768        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
    4869        fs(NComp)               as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
    49         T0                              as temperature  (Brief="Reference temperature", Default=298.15);
     70        To                              as temperature  (Brief="Reference temperature", Default=298.15);
    5071       
    5172        VARIABLES
    5273out Outlet                      as vapour_stream; # Outlet stream
    5374
    54         Gs(NComp)               as energy_mol   (Brief="Gibbs energy in standard state");
     75        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
    5576        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier");
     77        activ(NComp)    as Real                 (Brief="Activity", Lower=1e-20);
    5678        phi(NComp)              as fugacity             (Brief="Fugacity coefficient", Default=1);
    57         activ(NComp)    as Real                 (Brief="Activity", Lower=1e-20);
    58 
     79       
    5980        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
    6081        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Default=0);
     
    80101        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
    81102
    82         "Gibbs Energy of Formation"
    83         Gs = PP.IdealGasGibbsOfFormation(Outlet.T);
    84 
    85 #       "Gibbs Energy of Formation without Cp correction"
    86 #       Gs = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/T0+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/T0);
     103        "Gibbs free-energy of formation"
     104        G = PP.IdealGasGibbsOfFormation(Outlet.T);
     105
     106#       "Gibbs free-energy of formation without Cp correction"
     107#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);
    87108
    88109        for i in [1:NComp]
    89110        "Lagrangian multiplier"
    90                 Gs(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
     111                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
    91112
    92113          if (Outletm.z(i) > 0) then
     
    110131        activ = phi*Outlet.P*Outlet.z/fs;
    111132end
     133
     134
     135#*---------------------------------------------------------------------
     136*       only liquid phase
     137*--------------------------------------------------------------------*#
     138Model gibbs_liq as tank_liq
     139        ATTRIBUTES
     140        Pallete = true;
     141        Icon    = "icon/cstr";
     142        Brief   = "Model of a generic liquid-phase Gibbs CSTR";
     143        Info    = "
     144Requires the information of:
     145* number of elements
     146* matrix of elements (elements by compoments)
     147";
     148
     149        PARAMETERS
     150outer NElem                     as Integer              (Brief="Number of elements", Default=1);
     151        Rg                              as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
     152        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
     153        Ps                              as pressure             (Brief="Pressure of standard state", Default=1, DisplayUnit='atm');
     154        To                              as temperature  (Brief="Reference temperature", Default=298.15);
     155       
     156        VARIABLES
     157out Outlet                      as liquid_stream; # Outlet stream
     158
     159        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
     160        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier");
     161        activ(NComp)    as Real                 (Brief="Activity", Lower=0);
     162        gamma(NComp)    as fugacity             (Brief="Activity coefficient", Default=1);
     163
     164        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
     165        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Default=0);
     166        Fi(NComp)               as flow_mol             (Brief="Component molar flow rate");
     167
     168        EQUATIONS
     169        "Outlet stream"
     170        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
     171       
     172        "Mechanical equilibrium"
     173        Outlet.P = Outletm.P;
     174
     175        "Steady-state"
     176        Outlet.F = sum(Fi);
     177
     178        "Component molar flow rate"
     179        Fi = Outlet.F*Outlet.z;
     180       
     181        "Energy balance"
     182        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
     183
     184        "Element balance"
     185        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
     186       
     187        "Gibbs free-energy of formation"
     188        G = PP.IdealGasGibbsOfFormation(Outlet.T);
     189
     190#       "Gibbs free-energy of formation without Cp correction"
     191#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);
     192
     193        for i in [1:NComp]
     194        "Lagrangian multiplier"
     195                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
     196
     197          if (Outletm.z(i) > 0) then
     198            "Molar conversion"
     199            Fi(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
     200          else if (Outlet.z(i) > 0) then
     201                        "Molar conversion"
     202                                conv(i) = 1;    # ?
     203                        else
     204                        "Molar conversion"
     205                                conv(i) = 0;    # ?
     206                        end
     207          end
     208       
     209        end
     210       
     211        "Activity coefficient"
     212        gamma = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
     213       
     214        "Activity"
     215        activ = gamma*Outlet.z*exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
     216end
Note: See TracChangeset for help on using the changeset viewer.