Changeset 415 for trunk/eml/reactors

Timestamp:

Nov 26, 2007, 9:50:55 AM (16 years ago)

Author:

Rodolfo Rodrigues

Message:

Added liquid-phase calculation

File:

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

trunk/eml/reactors/equil.mso

@@ -17 +17 @@
 *----------------------------------------------------------------------
 *
-*
+*   Description:
+*       Thermodynamic equilibrium modeling of a reactor based on 
+*       equilibrium constants approach.
+*
+*   Assumptions:
+*       * thermodynamic equilibrium
+*               * steady-state
+*
+*       Specify:
+*               * inlet stream
+*               * stoichiometric matrix
+*               * equilibrium temperature
 *
 *----------------------------------------------------------------------
@@ -30 +41 @@
 *       only vapour phase
 *--------------------------------------------------------------------*#
-
 Model equil_vap as tank_vap
         ATTRIBUTES
         Pallete         = true;
         Icon            = "icon/cstr"; 
-        Brief           = "Model of a generic vapour-phase Equilibrium CSTR";
+        Brief           = "Model of a generic vapour-phase equilibrium CSTR";
         Info            = "
 Requires the information of:
 * number of reactions
-* matrix of stoichiometric coefficients (compounds, reaction)
+* matrix of stoichiometric coefficients (components by reactions)
 ";
 
@@ -47 +57 @@
         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');
-        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 of formation");
         K(NReac)        as Real                 (Brief="Equillibrium constant",Default=1.5);
         activ(NComp)as Real             (Brief="Activity",Default=0.2);
-        phi(NComp)      as fugacity             (Brief="Fugacity coefficient", Default=1);
         
         rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
@@ -72 +81 @@
         
         "Steady-state"
-        Outlet.F = Inlet.F+sum(sumt(stoic*extent));
-
-        "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);
-
+        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 = exp(-sumt(Gs*stoic)/(Rg*Outlet.T));
-        
-        for j in [1:NReac]
-        "Equilibrium"
                 K(j) = prod(activ^stoic(:,j));
         end
@@ -107 +117 @@
         end
 
-        "Fugacity coefficient"
-        phi = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z);
-
         "Activity"
-        activ = phi*Outlet.P*Outlet.z/fs;
+        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            = "
+Requires the information of:
+* number of reactions
+* matrix of stoichiometric coefficients (components by reactions)
+";
+
+        PARAMETERS
+        NReac           as Integer              (Brief="Number of reactions", Default=1);
+    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix");
+        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; # Outlet stream
+
+        G(NReac)        as enth_mol     (Brief="Gibbs free-energy of formation");
+        K(NReac)        as fraction             (Brief="Equillibrium constant");
+        activ(NComp)as Real             (Brief="Activity");
+        
+        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
+        extent(NReac)as flow_mol        (Brief="Extent of reaction");
+        conv(NComp) as Real             (Brief="Fractional conversion of component", Default=0);
+        
+        EQUATIONS
+        "Outlet stream"
+        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*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) > 0) 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