Ignore:
Timestamp:
Mar 20, 2009, 4:59:07 PM (14 years ago)
Author:
gerson bicca
Message:

PFR reactor updates

File:
1 edited

Legend:

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

    r574 r745  
    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                 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");
     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");
    4551
    4652SET
    4753       
    48         Mw   = PP.MolecularWeight();
     54        Mw              = PP.MolecularWeight();
     55        Across  = 0.25*pi*Dpipe^2;
     56        dx                      = L/NDisc;
     57        dv                      = Across*dx;
    4958
    5059VARIABLES
    5160
    52 in      Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
     61in              Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
    5362out     Outlet  as stream       (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}");
    5463
    55         str(NDisc+1) as streamPH;
    56         vel(NDisc+1) as velocity;
    57         vol(NDisc+1) as vol_mol;
    58         rho(NDisc+1) as dens_mass;
     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}");
    5971       
    60         q(NDisc)                                as heat_rate;
    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;
     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");
    6780
    6881EQUATIONS
     
    8194        Outlet.h = str(NDisc+1).h;
    8295        Outlet.v = str(NDisc+1).v;
     96
     97"Reactor Initial Length"
     98        Lincr(1) = 0*'m';
     99
     100        C(:,1)*vol(1) = Inlet.z;
    83101       
    84         for z in [1:NDisc] do
    85                 for c in [1:NComp-1] do
    86                         "Component Molar Balance"
    87                         diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
    88                                 + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
    89                 end
     102for i in [1:NDisc] do
    90103
    91                 "Energy Balance"
    92                 diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
    93                         sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);
     104for c in [1:NComp] do
    94105
    95                 "Energy Holdup"
    96                 E(z) = Mt(z)*str(z+1).h - str(z+1).P*Across*L/NDisc;
     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;
    97109
    98                 "mass flow is considered constant"
    99                 str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar)
    100                 #rho(z+1)*vel(z+1) = rho(z)*vel(z);  # FIXME: this is correct! But does not converge.
     110end
     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
     139end
     140
     141for 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;
    101160       
    102                 "Molar concentration"
    103                 C(:,z) * Across*L/NDisc = M(:,z);
    104                
    105                 "Sum of M"
    106                 Mt(z) = sum(M(:,z));
    107                
    108                 "Geometrical constraint"
    109                 Across*L/NDisc = Mt(z) * vol(z);
     161"Velocity"
     162        str(i).F = vel(i)*Across*rhom(i);
    110163
    111                 "Molar fraction"
    112                 str(z+1).z * Mt(z) = M(:,z);
    113         end
     164if 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;
    114174
    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);
     175end
    121176
    122                 "Velocity"
    123                 vel(z)*Across = str(z).F*vol(z);
    124         end
    125177end
     178
     179end
Note: See TracChangeset for help on using the changeset viewer.