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

Fixed PFR model and sample

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.