Changeset 813


Ignore:
Timestamp:
Aug 5, 2009, 4:00:16 PM (14 years ago)
Author:
Rafael de Pelegrini Soares
Message:

Fixed PFR model and sample

Location:
branches/gui
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/gui/eml/reactors/pfr.mso

    r745 r813  
    3535PARAMETERS
    3636
    37 outer   PP                                                      as Plugin               (Brief = "External Physical Properties", Type="PP");
     37outer PP                                                        as Plugin               (Brief = "External Physical Properties", Type="PP");
    3838outer   NComp                                           as Integer              (Brief="Number of components");
    39                
    40         NReac                                                   as Integer              (Brief="Number of reactions");
    41         stoic(NComp, NReac)     as Real                         (Brief = "Stoichiometric Matrix");
    42         NDisc                                                           as Integer              (Brief="Number of points of discretization", Default=10);
    43         Mw(NComp)                                       as molweight    (Brief="Component Mol Weight");
    44         L                                                                               as length                       (Brief="Reactor Length");
    45         Across                                                  as area                 (Brief="Cross section area");
    46         Dpipe                                                           as length                       (Brief="Reactor Inner Diameter");
    47         pi                                                                      as Real                         (Brief="pi number",Default=3.141592, Symbol = "\pi");
    48         dx                                                                      as length               (Brief = "Incremental Length");
    49         dv                                                                      as volume               (Brief = "Incremental volume");
    50         Roughness                                       as length                       (Brief="Reactor Tube Roughness", Default = 4.572E-5, Symbol = "\varepsilon");
     39                NReac                                           as Integer              (Brief="Number of reactions");
     40                stoic(NComp, NReac)     as Real                         (Brief = "Stoichiometric Matrix");
     41                NDisc                                           as Integer              (Brief="Number of points of discretization", Default=10);
     42                Mw(NComp)                               as molweight    (Brief="Component Mol Weight");
     43                L                                                               as length               (Brief="Reactor Length");
     44                Across                                          as area                 (Brief="Cross section area");
    5145
    5246SET
    5347       
    54         Mw              = PP.MolecularWeight();
    55         Across  = 0.25*pi*Dpipe^2;
    56         dx                      = L/NDisc;
    57         dv                      = Across*dx;
     48        Mw   = PP.MolecularWeight();
    5849
    5950VARIABLES
    6051
    61 in              Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
     52in      Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
    6253out     Outlet  as stream       (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}");
    6354
    64         str(NDisc+1)                    as streamPH;
    65         vel(NDisc+1)                    as velocity;
    66         vol(NDisc+1)                    as vol_mol;
    67         rho(NDisc+1)                    as dens_mass;
    68         rhom(NDisc+1)           as dens_mol;
    69         dP(NDisc+1)                     as press_delta          (Brief = "Friction Pressure Drop");
    70         Lincr(NDisc+1)          as length                       (Brief = "Length Points", Symbol = "L_{incr}");
     55        str(NDisc+1) as vapour_stream;
     56        vol(NDisc+1) as vol_mol;
     57        rho(NDisc+1) as dens_mass;
    7158       
    72         q(NDisc)                                                as heat_rate;
    73         C(NComp, NDisc+1)       as conc_mol                     (Brief="Components concentration");
    74         E(NDisc)                                                as energy                               (Brief="Total Energy Holdup on element");
    75         r(NReac, NDisc)                         as reaction_mol;
    76         Hr(NReac, NDisc)                        as heat_reaction;
    77         Re(NDisc+1)                             as Real                                 (Brief = "Reynolds Number Profile",Lower = 1E-6);
    78         mu(NDisc+1)                             as viscosity                    (Brief = "Viscosity Profile" , Symbol = "\mu");
    79         fns(NDisc+1)                                    as fricfactor           (Brief = "No Slip Friction Factor");
     59        q(NDisc)                                as heat_rate;
     60
     61        M(NComp, NDisc)         as mol                  (Brief = "Molar holdup");
     62        Mt(NDisc)                       as mol                  (Brief = "Molar holdup");
     63        C(NComp, NDisc)         as conc_mol     (Brief="Components concentration", Lower=-1e-6);
     64        E(NDisc)                                as energy               (Brief="Total Energy Holdup on element");
     65        r(NReac, NDisc)                 as reaction_mol;
     66        Hr(NReac, NDisc)        as heat_reaction;
     67        #Hf(NComp, NDisc)       as heat_reaction;
    8068
    8169EQUATIONS
     
    9482        Outlet.h = str(NDisc+1).h;
    9583        Outlet.v = str(NDisc+1).v;
     84       
     85        for z in [1:NDisc] do
     86                for c in [1:NComp] do
     87                        "Component Molar Balance"
     88                        diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
     89                                + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
     90                end
     91               
     92                "Energy Balance"
     93                diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
     94                        sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);
    9695
    97 "Reactor Initial Length"
    98         Lincr(1) = 0*'m';
     96                "Energy Holdup"
     97                E(z) = Mt(z)*str(z+1).h - str(z+1).P*Across*L/NDisc;
    9998
    100         C(:,1)*vol(1) = Inlet.z;
     99                "mass flow is considered constant"
     100                str(z+1).F*vol(z+1)*rho(z+1) = str(z).F*vol(z)*rho(z);
    101101       
    102 for i in [1:NDisc] do
     102                "Molar concentration"
     103                C(:,z) * Across*L/NDisc = M(:,z);
     104               
     105                "Geometrical constraint"
     106                Across*L/NDisc = Mt(z) * vol(z);
    103107
    104 for c in [1:NComp] do
     108                "Molar fraction"
     109                str(z+1).z * Mt(z) = M(:,z);
     110               
     111                #Hf(:,z) = PP.IdealGasEnthalpyOfFormation(str(z+1).T);
     112                #Hr(:,z) = -sum(stoic*Hf(:, z));
     113        end
    105114
    106 "Component Molar Balance"
    107         #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;
    108         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;
    109 
     115        for z in [1:NDisc+1] do
     116                "Specific Volume"
     117                vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z);
     118               
     119                "Specific Mass"
     120                rho(z) = PP.VapourDensity(str(z).T, str(z).P, str(z).z);
     121        end
    110122end
    111 
    112 "Constraint"
    113         sum(str(i+1).z ) = 1;
    114 #       Mt(i) = sum(M(:,i));
    115 
    116 "Molar concentration"
    117         #C(:,i) = rhom(i)*str(i).z;
    118         str(i+1).z = C(:,i+1) * vol(i+1);
    119 
    120 #"Molar fraction"
    121         #str(i+1).z * Mt(i) = M(:,i);
    122 
    123 #"Geometrical constraint"
    124         #Mt(i)  = dv * rhom(i);
    125 
    126 "Energy Balance"
    127         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);
    128 
    129 "Energy Holdup"
    130         #E(i) = Mt(i)*str(i+1).h - str(i+1).P*dv;
    131         E(i) = dv/vol(i+1) * str(i+1).h - str(i+1).P*dv;
    132 
    133 "Outlet Pressure"
    134         str(i+1).P =str(1).P - dP(i+1);
    135 
    136 "Incremental Length"
    137         Lincr(i+1) = i*dx;
    138 
    139 end
    140 
    141 for i in [1:NDisc+1] do
    142 
    143 "Incremental Pressure Drop"
    144         dP(i) = 0.5*fns(i)*Lincr(i)*rho(i)*vel(i)*vel(i)/Dpipe; # simple equation for pressure drop (Darcy Equation)
    145 
    146 "Specific Volume"
    147         vol(i) = PP.VapourVolume(str(i).T, str(i).P, str(i).z);
    148 
    149 "Specific Mass"
    150         rho(i) = PP.VapourDensity(str(i).T, str(i).P, str(i).z);
    151 
    152 "Molar Density"
    153         rhom(i)* vol(i) = 1;
    154 
    155 "Viscosity"
    156         mu(i) = 0.027*'cP';#PP.VapourViscosity(str(i).T,str(i).P,str(i).z); # VRTherm mu=16 cP !!!!!!
    157 
    158 "Reynolds Number"
    159         Re(i)*mu(i)     = rho(i)*vel(i)*Dpipe;
    160        
    161 "Velocity"
    162         str(i).F = vel(i)*Across*rhom(i);
    163 
    164 if Re(i) > 2300
    165        
    166         then
    167 "Friction Factor  for Pressure Drop - Turbulent Flow"
    168         #1/sqrt(fns(i))= -2*log(abs(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i))));
    169         1/sqrt(fns(i))= -2*log(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i)));
    170        
    171         else
    172 "Friction Factor for Pressure Drop - laminar Flow"
    173         fns(i)*Re(i) = 16;
    174 
    175 end
    176 
    177 end
    178 
    179 end
  • branches/gui/sample/reactors/sample_pfr.mso

    r745 r813  
    3939PARAMETERS
    4040
    41         PP              as Plugin (Brief="Physical Properties",
     41        PP as Plugin (Brief="Physical Properties",
    4242                Type = "PP",
    43                 Components = ["acetone", "acetic anhydride", "methane" ],
     43                Components = ["acetone", "ketene", "methane" ],
    4444                LiquidModel = "PR",
    4545                VapourModel = "PR"
    4646        );
    4747        NComp   as Integer;
     48SET
     49        NComp = PP.NumberOfComponents;
    4850
    4951DEVICES
     52        s1      as source;
     53        Reac    as pfr;
    5054
    51         s1      as simple_source;
    52         Reac as pfr;
     55CONNECTIONS
     56        s1.Outlet  to  Reac.Inlet;
     57       
     58SET
     59        Reac.NDisc      = 15;
     60        Reac.L          = 2.28 * 'm';
     61        # Area was chosen to match the original example data
     62        Reac.Across = 1.27*'m^3'/Reac.L;
    5363
    54 SET
    55         NComp = PP.NumberOfComponents;
    56        
    57         s1.ValidPhases = "Vapour-Liquid";
    58         #s1.CompositionBasis = "Molar";
    59        
    60         Reac.NDisc      = 10;
    61         Reac.Dpipe      = 3 * 'in';
    62         Reac.L                  = 15 * 'm';
     64        # Only one reaction
    6365        Reac.NReac      = 1;
    64         Reac.Roughness  = 4.572E-5*'m';
     66        # Reaction 1: acetone -> ketene + methane
     67        Reac.stoic(:,1) = [-1, 1, 1];
    6568
    66         # Reaction 1: A -> B + C
    67         Reac.stoic(:,1) = [-1, 1, 1];
     69SPECIFY
     70        # Feed specification
     71        s1.Fw = 8000 * 'kg/h';
     72        s1.T = 1035 * 'K';
     73        s1.P = 1.6  * 'atm';
     74        s1.Composition = [0.98, 0.01, 0.01];
     75
     76        # Adiabatic
     77        Reac.q     = 0     * 'J/s';
    6878       
    6979EQUATIONS
     
    7787        Reac.Hr(1,z) = -80.77 * 'kJ/mol';
    7888
     89"Pressure Drop (no pressure drop)"
     90        Reac.str(z+1).P = Reac.str(z).P;
     91
    7992end
    8093
    81 SPECIFY
    82         s1.F = 138/1000 * 'kmol/h';
    83         s1.T = 1035 * 'K';
    84         s1.P = 1.6  * 'atm';
    85         s1.MolarComposition = [1, 0, 0];
    86 
    87         Reac.q     = 0     * 'J/s';
    88        
    89 CONNECTIONS
    90 
    91         s1.Outlet  to  Reac.Inlet;
    92        
    9394INITIAL
    9495
    9596for z in [2:Reac.NDisc+1] do
     97        Reac.str(z).T = Reac.Inlet.T;
     98        Reac.str(z).z(1:NComp) = Reac.Inlet.z(1:NComp);
     99end
    96100       
    97         Reac.str(z).T = Reac.Inlet.T;
    98         Reac.vel(z) = Reac.vel(1);
    99         #Reac.str(z).F = Reac.str(1).F;
    100         Reac.str(z).z(1:NComp-1) = Reac.Inlet.z(1:NComp-1);
    101 
    102 end
    103 
    104 VARIABLES
    105 
    106 SOMA as Real;
    107 
    108 EQUATIONS
    109 
    110 SOMA = sum(Reac.Outlet.z);
    111 
    112101OPTIONS
    113         # This model can be solved in both steady or dynamic
    114102        Dynamic         = true;
    115103        TimeStep        = 0.05;
    116         TimeEnd         = 3;
    117         DAESolver(File="sundials", RelativeAccuracy=1e-2);
    118         NLASolver(RelativeAccuracy=1e-4);
     104        TimeEnd         = 10;
    119105        #InitialFile = "PFR_AceticAnhydride";
    120         #Dynamic = false;
    121         #DAESolver = "dassl";
     106        #Dynamic = false; # Requires a GuessFile from a dynamic simulation
    122107end
Note: See TracChangeset for help on using the changeset viewer.