Changeset 745

Timestamp:

Mar 20, 2009, 4:59:07 PM (13 years ago)

Author:

gerson bicca

Message:

PFR reactor updates

Location:

branches/gui

Files:

• eml/reactors/pfr.mso (modified) (2 diffs)
• sample/reactors/PFR_AceticAnhydride.rlt (added)
• sample/reactors/sample_pfr.mso (modified) (6 diffs)

branches/gui/eml/reactors/pfr.mso

@@ -35 +35 @@
 PARAMETERS
 
-outer PP                                                        as Plugin               (Brief = "External Physical Properties", Type="PP");
+outer   PP                                                      as Plugin               (Brief = "External Physical Properties", Type="PP");
 outer   NComp                                           as Integer              (Brief="Number of components");
-                NReac                                           as Integer              (Brief="Number of reactions");
-                stoic(NComp, NReac)     as Real                         (Brief = "Stoichiometric Matrix");
-                NDisc                                           as Integer              (Brief="Number of points of discretization", Default=10);
-                Mw(NComp)                               as molweight    (Brief="Component Mol Weight");
-                L                                                               as length               (Brief="Reactor Length");
-                Across                                          as area                 (Brief="Cross section area");
+                
+        NReac                                                   as Integer              (Brief="Number of reactions");
+        stoic(NComp, NReac)     as Real                         (Brief = "Stoichiometric Matrix");
+        NDisc                                                           as Integer              (Brief="Number of points of discretization", Default=10);
+        Mw(NComp)                                       as molweight    (Brief="Component Mol Weight");
+        L                                                                               as length                       (Brief="Reactor Length");
+        Across                                                  as area                 (Brief="Cross section area");
+        Dpipe                                                           as length                       (Brief="Reactor Inner Diameter");
+        pi                                                                      as Real                         (Brief="pi number",Default=3.141592, Symbol = "\pi");
+        dx                                                                      as length               (Brief = "Incremental Length");
+        dv                                                                      as volume               (Brief = "Incremental volume");
+        Roughness                                       as length                       (Brief="Reactor Tube Roughness", Default = 4.572E-5, Symbol = "\varepsilon");
 
 SET
         
-        Mw   = PP.MolecularWeight();
+        Mw              = PP.MolecularWeight();
+        Across  = 0.25*pi*Dpipe^2; 
+        dx                      = L/NDisc;
+        dv                      = Across*dx;
 
 VARIABLES
 
-in      Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
+in              Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
 out     Outlet  as stream       (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}");
 
-        str(NDisc+1) as streamPH;
-        vel(NDisc+1) as velocity;
-        vol(NDisc+1) as vol_mol;
-        rho(NDisc+1) as dens_mass;
+        str(NDisc+1)                    as streamPH;
+        vel(NDisc+1)                    as velocity;
+        vol(NDisc+1)                    as vol_mol;
+        rho(NDisc+1)                    as dens_mass;
+        rhom(NDisc+1)           as dens_mol;
+        dP(NDisc+1)                     as press_delta          (Brief = "Friction Pressure Drop");
+        Lincr(NDisc+1)          as length                       (Brief = "Length Points", Symbol = "L_{incr}");
         
-        q(NDisc)                                as heat_rate;
-        M(NComp, NDisc)         as mol                  (Brief = "Molar holdup");
-        Mt(NDisc)                       as mol                  (Brief = "Molar holdup");
-        C(NComp, NDisc)         as conc_mol     (Brief="Components concentration", Lower=-1e-6);
-        E(NDisc)                                as energy               (Brief="Total Energy Holdup on element");
-        r(NReac, NDisc)                 as reaction_mol;
-        Hr(NReac, NDisc)        as heat_reaction;
+        q(NDisc)                                                as heat_rate;
+        C(NComp, NDisc+1)       as conc_mol                     (Brief="Components concentration");
+        E(NDisc)                                                as energy                               (Brief="Total Energy Holdup on element");
+        r(NReac, NDisc)                         as reaction_mol;
+        Hr(NReac, NDisc)                        as heat_reaction;
+        Re(NDisc+1)                             as Real                                 (Brief = "Reynolds Number Profile",Lower = 1E-6);
+        mu(NDisc+1)                             as viscosity                    (Brief = "Viscosity Profile" , Symbol = "\mu");
+        fns(NDisc+1)                                    as fricfactor           (Brief = "No Slip Friction Factor");
 
 EQUATIONS
@@ -81 +94 @@
         Outlet.h = str(NDisc+1).h;
         Outlet.v = str(NDisc+1).v;
+
+"Reactor Initial Length"
+        Lincr(1) = 0*'m';
+
+        C(:,1)*vol(1) = Inlet.z;
         
-        for z in [1:NDisc] do
-                for c in [1:NComp-1] do
-                        "Component Molar Balance"
-                        diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
-                                + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
-                end
+for i in [1:NDisc] do
 
-                "Energy Balance"
-                diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
-                        sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);
+for c in [1:NComp] do
 
-                "Energy Holdup"
-                E(z) = Mt(z)*str(z+1).h - str(z+1).P*Across*L/NDisc;
+"Component Molar Balance"
+        #diff(M(c,i)) = str(i).F*str(i).z(c) - str(i+1).F*str(i+1).z(c) + sum(stoic(c,:)*r(:, i)) * dv;
+        diff(C(c,i+1)*dv) = (vel(i)*C(c,i) - vel(i+1)*C(c,i+1) )*Across + sum(stoic(c,:)*r(:, i))*dv;
 
-                "mass flow is considered constant"
-                str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar)
-                #rho(z+1)*vel(z+1) = rho(z)*vel(z);  # FIXME: this is correct! But does not converge.
+end
+
+"Constraint"
+        sum(str(i+1).z ) = 1;
+#       Mt(i) = sum(M(:,i));
+
+"Molar concentration"
+        #C(:,i) = rhom(i)*str(i).z;
+        str(i+1).z = C(:,i+1) * vol(i+1);
+
+#"Molar fraction"
+        #str(i+1).z * Mt(i) = M(:,i);
+
+#"Geometrical constraint"
+        #Mt(i)  = dv * rhom(i);
+
+"Energy Balance"
+        diff(E(i)) = str(i).F*str(i).h - str(i+1).F*str(i+1).h +        sum(Hr(:,i)*r(:,i)) * dv - q(i);
+
+"Energy Holdup"
+        #E(i) = Mt(i)*str(i+1).h - str(i+1).P*dv;
+        E(i) = dv/vol(i+1) * str(i+1).h - str(i+1).P*dv;
+
+"Outlet Pressure"
+        str(i+1).P =str(1).P - dP(i+1);
+
+"Incremental Length"
+        Lincr(i+1) = i*dx;
+
+end
+
+for i in [1:NDisc+1] do
+
+"Incremental Pressure Drop"
+        dP(i) = 0.5*fns(i)*Lincr(i)*rho(i)*vel(i)*vel(i)/Dpipe; # simple equation for pressure drop (Darcy Equation)
+
+"Specific Volume"
+        vol(i) = PP.VapourVolume(str(i).T, str(i).P, str(i).z);
+
+"Specific Mass"
+        rho(i) = PP.VapourDensity(str(i).T, str(i).P, str(i).z);
+
+"Molar Density"
+        rhom(i)* vol(i) = 1;
+
+"Viscosity" 
+        mu(i) = 0.027*'cP';#PP.VapourViscosity(str(i).T,str(i).P,str(i).z); # VRTherm mu=16 cP !!!!!!
+
+"Reynolds Number"
+        Re(i)*mu(i)     = rho(i)*vel(i)*Dpipe;
         
-                "Molar concentration"
-                C(:,z) * Across*L/NDisc = M(:,z);
-                
-                "Sum of M"
-                Mt(z) = sum(M(:,z));
-                
-                "Geometrical constraint"
-                Across*L/NDisc = Mt(z) * vol(z);
+"Velocity"
+        str(i).F = vel(i)*Across*rhom(i);
 
-                "Molar fraction"
-                str(z+1).z * Mt(z) = M(:,z);
-        end
+if Re(i) > 2300
+        
+        then
+"Friction Factor  for Pressure Drop - Turbulent Flow"
+        #1/sqrt(fns(i))= -2*log(abs(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i))));
+        1/sqrt(fns(i))= -2*log(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i)));
+        
+        else
+"Friction Factor for Pressure Drop - laminar Flow"
+        fns(i)*Re(i) = 16;
 
-        for z in [1:NDisc+1] do
-                "Specific Volume"
-                vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z);
-                
-                "Specific Mass"
-                rho(z) = PP.VapourDensity(str(z).T, str(z).P, str(z).z);
+end
 
-                "Velocity"
-                vel(z)*Across = str(z).F*vol(z);
-        end
 end
+
+end

branches/gui/sample/reactors/sample_pfr.mso

@@ -49 +49 @@
 DEVICES
 
-        s1      as source;
+        s1      as simple_source;
         Reac as pfr;
 
@@ -56 +56 @@
         
         s1.ValidPhases = "Vapour-Liquid";
-        s1.CompositionBasis = "Molar";
+        #s1.CompositionBasis = "Molar";
         
         Reac.NDisc      = 10;
-        Reac.Across     = 0.7 * 'in^2';
-        Reac.L                  = 2.28 * 'm';
+        Reac.Dpipe      = 3 * 'in';
+        Reac.L                  = 15 * 'm';
         Reac.NReac      = 1;
+        Reac.Roughness  = 4.572E-5*'m';
 
         # Reaction 1: A -> B + C
@@ -76 +77 @@
         Reac.Hr(1,z) = -80.77 * 'kJ/mol';
 
-"Pressure Drop (no pressure drop)"
-        Reac.str(z+1).P = Reac.str(z).P;
-
 end
 
@@ -85 +83 @@
         s1.T = 1035 * 'K';
         s1.P = 1.6  * 'atm';
-        s1.Composition = [1, 0, 0];
+        s1.MolarComposition = [1, 0, 0];
 
         Reac.q     = 0     * 'J/s';
@@ -98 +96 @@
         
         Reac.str(z).T = Reac.Inlet.T;
+        Reac.vel(z) = Reac.vel(1);
+        #Reac.str(z).F = Reac.str(1).F;
         Reac.str(z).z(1:NComp-1) = Reac.Inlet.z(1:NComp-1);
 
 end
-        
+
+VARIABLES
+
+SOMA as Real;
+
+EQUATIONS
+
+SOMA = sum(Reac.Outlet.z);
+
 OPTIONS
         # This model can be solved in both steady or dynamic
@@ -107 +115 @@
         TimeStep        = 0.05;
         TimeEnd         = 3;
+        DAESolver(File="sundials", RelativeAccuracy=1e-2);
+        NLASolver(RelativeAccuracy=1e-4);
+        #InitialFile = "PFR_AceticAnhydride";
         #Dynamic = false;
         #DAESolver = "dassl";