Changeset 414 for trunk/eml/reactors

Timestamp:

Nov 25, 2007, 2:15:43 PM (16 years ago)

Author:

Rodolfo Rodrigues

Message:

Added liquid-phase calculation

File:

• trunk/eml/reactors/gibbs.mso (modified) (5 diffs)

trunk/eml/reactors/gibbs.mso

@@ -17 +17 @@
 *----------------------------------------------------------------------
 *
-*
+*   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
 *
 *----------------------------------------------------------------------
@@ -30 +52 @@
 *       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            = "
+        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, compounds)
+* matrix of elements (elements by compoments)
 ";
 
@@ -47 +68 @@
         na(NElem,NComp) as Real                 (Brief="Number of elements per component");
         fs(NComp)               as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
-        T0                              as temperature  (Brief="Reference temperature", Default=298.15);
+        To                              as temperature  (Brief="Reference temperature", Default=298.15);
         
         VARIABLES
 out Outlet                      as vapour_stream; # Outlet stream
 
-        Gs(NComp)               as energy_mol   (Brief="Gibbs energy in standard state");
+        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);
-        activ(NComp)    as Real                 (Brief="Activity", Lower=1e-20);
-
+        
         rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
         conv(NComp)     as Real                 (Brief="Fractional conversion of component", Default=0);
@@ -80 +101 @@
         sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
 
-        "Gibbs Energy of Formation"
-        Gs = PP.IdealGasGibbsOfFormation(Outlet.T);
-
-#       "Gibbs Energy of Formation without Cp correction"
-#       Gs = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/T0+PP.IdealGasEnthalpyOfFormationAt25C()*(1-Outlet.T/T0);
+        "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"
-                Gs(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
+                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
 
           if (Outletm.z(i) > 0) then
@@ -110 +131 @@
         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