Changeset 415 for trunk/eml


Ignore:
Timestamp:
Nov 26, 2007, 9:50:55 AM (15 years ago)
Author:
Rodolfo Rodrigues
Message:

Added liquid-phase calculation

File:
1 edited

Legend:

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

    r403 r415  
    1717*----------------------------------------------------------------------
    1818*
    19 *
     19*   Description:
     20*       Thermodynamic equilibrium modeling of a reactor based on
     21*       equilibrium constants approach.
     22*
     23*   Assumptions:
     24*       * thermodynamic equilibrium
     25*               * steady-state
     26*
     27*       Specify:
     28*               * inlet stream
     29*               * stoichiometric matrix
     30*               * equilibrium temperature
    2031*
    2132*----------------------------------------------------------------------
     
    3041*       only vapour phase
    3142*--------------------------------------------------------------------*#
    32 
    3343Model equil_vap as tank_vap
    3444        ATTRIBUTES
    3545        Pallete         = true;
    3646        Icon            = "icon/cstr";
    37         Brief           = "Model of a generic vapour-phase Equilibrium CSTR";
     47        Brief           = "Model of a generic vapour-phase equilibrium CSTR";
    3848        Info            = "
    3949Requires the information of:
    4050* number of reactions
    41 * matrix of stoichiometric coefficients (compounds, reaction)
     51* matrix of stoichiometric coefficients (components by reactions)
    4252";
    4353
     
    4757        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
    4858        fs(NComp)       as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
    49         T0                              as temperature  (Brief="Reference temperature", Default=298.15);
     59        To                      as temperature  (Brief="Reference temperature", Default=298.15);
    5060
    5161        VARIABLES
    5262out Outlet              as vapour_stream; # Outlet stream
    5363
    54         Gs(NComp)       as energy_mol   (Brief="Gibbs energy in standard state");
     64        G(NComp)        as energy_mol   (Brief="Gibbs free-energy of formation");
    5565        K(NReac)        as Real                 (Brief="Equillibrium constant",Default=1.5);
    5666        activ(NComp)as Real             (Brief="Activity",Default=0.2);
    57         phi(NComp)      as fugacity             (Brief="Fugacity coefficient", Default=1);
    5867       
    5968        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
     
    7281       
    7382        "Steady-state"
    74         Outlet.F = Inlet.F+sum(sumt(stoic*extent));
    75 
    76         "Gibbs Energy of Formation"
    77         Gs = PP.IdealGasGibbsOfFormation(Outlet.T);
    78 
    79 #       "Gibbs Energy of Formation without Cp correction"
    80 #       Gs = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/T0+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/T0);
    81 
     83        Outlet.F = Inlet.F + sum(sumt(stoic*extent));
     84
     85        "Gibbs free-energy of formation"
     86        G = PP.IdealGasGibbsOfFormation(Outlet.T);
     87
     88#       "Gibbs free-energy of formation without Cp correction"
     89#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);
     90
     91        "Gibbs free energy of reaction"
     92#       sumt(G*stoic) = -Rg*Outlet.T*ln(K);
     93        K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
     94       
     95        for j in [1:NReac]
    8296        "Equilibrium constant"
    83         K = exp(-sumt(Gs*stoic)/(Rg*Outlet.T));
    84        
    85         for j in [1:NReac]
    86         "Equilibrium"
    8797                K(j) = prod(activ^stoic(:,j));
    8898        end
     
    107117        end
    108118
    109         "Fugacity coefficient"
    110         phi = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
    111 
    112119        "Activity"
    113         activ = phi*Outlet.P*Outlet.z/fs;
     120        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
    114121end
     122
     123
     124#*---------------------------------------------------------------------
     125*       only liquid-phase
     126*--------------------------------------------------------------------*#
     127Model equil_liq as tank_liq
     128        ATTRIBUTES
     129        Pallete         = true;
     130        Icon            = "icon/cstr";
     131        Brief           = "Model of a generic liquid-phase equilibrium CSTR";
     132        Info            = "
     133Requires the information of:
     134* number of reactions
     135* matrix of stoichiometric coefficients (components by reactions)
     136";
     137
     138        PARAMETERS
     139        NReac           as Integer              (Brief="Number of reactions", Default=1);
     140    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix");
     141        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
     142        Ps                      as pressure     (Brief="Standard pressure", Default=1, DisplayUnit='bar');
     143        To                      as temperature  (Brief="Reference temperature", Default=298.15);
     144
     145        VARIABLES
     146out Outlet              as liquid_stream; # Outlet stream
     147
     148        G(NReac)        as enth_mol     (Brief="Gibbs free-energy of formation");
     149        K(NReac)        as fraction             (Brief="Equillibrium constant");
     150        activ(NComp)as Real             (Brief="Activity");
     151       
     152        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
     153        extent(NReac)as flow_mol        (Brief="Extent of reaction");
     154        conv(NComp) as Real             (Brief="Fractional conversion of component", Default=0);
     155       
     156        EQUATIONS
     157        "Outlet stream"
     158        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V;
     159       
     160        "Mechanical equilibrium"
     161        Outlet.P = Outletm.P;
     162       
     163        "Energy balance"
     164        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
     165       
     166        "Steady-state"
     167        Outlet.F = Inlet.F + sum(sumt(stoic*extent));
     168       
     169        "Gibbs free-energy of formation"
     170        G = PP.IdealGasGibbsOfFormation(Outlet.T);
     171
     172#       "Gibbs free-energy of formation without Cp correction"
     173#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/To);
     174
     175        "Gibbs free energy of reaction"
     176#       sumt(G*stoic) = -Rg*Outlet.T*ln(K);
     177        K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
     178       
     179        for j in [1:NReac]
     180        "Equilibrium constant"
     181                K(j) = prod(activ^stoic(:,j));
     182        end
     183       
     184        for i in [1:NComp]
     185        "Outlet molar fraction"
     186                Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
     187        end     
     188
     189        for i in [1:NComp]
     190          if (Outletm.z(i) > 0) then
     191            "Molar conversion"
     192            Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
     193          else if (Outlet.z(i) > 0) then
     194                        "Molar conversion"
     195                                conv(i) = 1;    # ?
     196                        else
     197                        "Molar conversion"
     198                                conv(i) = 0;    # ?
     199                        end
     200          end
     201        end
     202       
     203        "Activity"
     204        activ = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z
     205                *exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
     206end
Note: See TracChangeset for help on using the changeset viewer.